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SUMMARY 


A modern control theory design package (ORACLS) for constructing control- 
lers and optimal filters for systems modeled by linear time-invariant differen- 
tial or difference equations is described. The digital FORTRAN-coded ORACLS 
system represents an application of some of today’s best numerical linear- 
algebra procedures to implement the linear-quadratic-Gaussian (LQG) methodology 
of modern control theory. Included are algorithms for computing eigensystems 
of real matrices, the relative stability of a matrix, factored forms for non- 
negative definite matrices, the solutions and least squares approximations to 
the solutions of certain linear matrix algebraic equations, the controllability 
properties of a linear time-invariant system, and the steady-state covariance 
matrix of an open-loop stable system forced by white noise. These subroutines 
are applied to implement the various techniques of the LQG methodology. Sub- 
routines are provided for solving both the continuous and discrete optimal linear 
regulator problems with noise-free measurements and the sampled-data optimal 
linear regulator problem. For measurement noise, duality theory and the optimal 
regulator algorithms are used to solve the continuous and discrete Kalman-Bucy 
filter problems. Subroutines are also included which give control laws causing 
the output of a system to track the output of a prescribed model. Finally, 
numerical examples are presented to illustrate the capability of the ORACLS 
system. 


INTRODUCTION 

Optimal linear quadratic regulator theory, currently referred to as the 
linear-quadratic-Gaussian (LQG) problem (ref. 1), has become the most widely 
accepted method of determining optimal control policy. Within this class of 
optimal control problems, the infinite-duration time-invariant version, which 
leads to constant-gain feedback control laws and constant Kalman-Bucy filter 
gains for reconstruction of the system state, has received the bulk of the 
attention because of its tractability and potential ease of implementation. 

The theory is particularly attractive to the control system designer because 
it provides a rigorous tool for dealing with multi-input and multi-output 
dynamical systems in both continuous and discrete form. 

During the same period of time that the LQG methodology was being devel- 
oped, there also appeared in the field of numerical linear algebra a variety 
of new and more efficient methods for analyzing the types of equations which 
occur in the time-invariant formulation of the LQG problem (ref. 2). The best 
of these numerical methods (in the opinion of the author) have been incorpo- 
rated into a modern efficient digital computer package which provides solutions 
to time- invariant continuous or discrete LQG problems with equal ease. The 
package is entitled Optimal Regulator Algorithms for the Control of Linear 
Systems (ORACLS). 



An overview of the features of the ORACLS system is presented in the next 
section. Following sections give a detailed description of the package con- 
tents along with numerical examples to illustrate the use of ORACLS to solve 
selected LQG problems. Software for the ORACLS system may be obtained from the 
centralized facility COSMIC located at the Computer Software Management and 
Information Center, Suite 112, Barrow Hall, University of Georgia, Athens, 

GA 30602, by requesting program LAR-12313. 


OVERVIEW OF ORACLS 

The ORACLS programing system is a collection of FORTRAN-coded subroutines 
to formulate, manipulate, and solve various LQG design problems. In order to 
apply ORACLS, the user is required to provide an executive (driver) program 
which inputs the problem coefficients, formulates and selects the system sub- 
routines to be used to solve a particular optimal control problem, and outputs 
desired information. ORACLS is constructed to allow the user considerable 
flexibility at each operational state. This flexibility is accomplished by pro- 
viding primary -subroutines at four levels: input-output, basic vector-matrix 

operations, analysis of linear time-invariant systems, and control synthesis 
based on LQG methodology. ORACLS provides a means of controlling program size 
by employing dynamic (vector) data storage. Except for certain subroutines 
obtained from other sources, data arrays in all ORACLS subroutines are treated 
as packed one-dimensional arrays which can easily be passed between subroutines 
without a maximum array size parameter appearing as an argument of the calling 
sequence. This dynamic storage capability allows program size to be specified 
and controlled through the user's driver program. In addition, ORACLS only 
loads those programs from the library which are called by the executive pro- 
gram, making the total machine requirements very flexible. As a result, ORACLS 
can be made to execute efficiently on a wide variety of computing machinery. 

The input-output category of ORACLS has subroutines for inputing (subrou- 
tine READ) and outputing (PRNT) numerical matrices. Additional subroutines 
allow for printing header information (RDTITL) and accumulation of output line- 
count information (LNCNT). 

The next category has subroutines for the basic vector-matrix operations 
of equation (EQUATE), addition (ADD), subtraction (SUBT), and multiplication 
(MULT). It also contains routines for scaling (SCALE), juxtaposition (JUXTC 
and JUXTR), and construction of matrix norms (MAXEL and NORMS), trace (TRCE), 
transpose (TRANP), and null and identity matrices (NULL and UNITY). 

The analysis category provides special and general purpose algorithms for 
computing (1) eigenvalues and eigenvectors of real matrices (EIGEN) by using 
the QR algorithm (ref. 2), (2) the relative stability of a given matrix 
(TESTSTA) , (3) matrix factorization (FACTOR), (4) the solution of linear 
constant-coefficient vector-matrix algebraic equations (SYMPDS, GELIM, and 
SNVDEC), (5) the controllability properties of a linear time-invariant system 
(CTROL), (6) the steady-state covariance matrix of an open-loop stable system 
forced by white noise (VARANCE), and (7) the transient response of continuous 
linear time- invariant systems (TRANSIT). 
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Algorithms are provided for the solution of real matrix equations. In the 
equation, 


AX = B 


( 1 ) 


with A positive definite, the Cholesky decomposition method (ref. 2) is 
applied (SYMPDS). Gaussian elimination (LU factorization from ref. 2) is used 
for the general case of nonsingular A (GELIM)- For rectangular or singular 
matrices A, singular-value decomposition procedures found in subroutine SNVDEC 
(ref. 2) can be used to find a solution of equation (1) in the least squares 
sense, to compute the pseudoinverse of A, and to find an orthonormal basis for 
the range space of A. A maximal rank matrix factorization for a positive 
semidefinite matrix can be obtained from FACTOR. 

For solution of the matrix equation, 


X = AXB + C 


( 2 ) 


the contraction mapping principle (ref. 3) is applied when A and B are 
asymptotically stable in the discrete sense (SUM). Equation (2) is used when 
solving the discrete steady-state Riccati equation by Newton’s method and also 
in the computation of the steady-state covariance matrix for a linear asymptoti- 
cally stable discrete system forced by white noise (ref. 4). 

Two subroutines are included for solving the matrix equation, 


AX + XB = C (3) 


For A and B admitting a unique solution X, the method of Bartels and 
Stewart (ref. 5) is used (BARSTW). For A and B asymptotically stable in 
the continuous sense, a subroutine implementing the bilinear transformation 
approach (ref. 6) is also included (BILIN). Equation (3) is used in solving 
the steady-state continuous matrix Riccati equation by Newton’s method, in 
finding the covariance statistics for continuous time-invariant systems forced 
by white noise, and in gain computation for observer theory (ref. 7). 

A subroutine (CTROL) is provided which computes the controllability matrix 
for a linear time-invariant dynamical system and, if this matrix is found to 
be rank deficient, also computes the system’s controllability canonical form 
(ref. 4) by application of the singular-value decomposition algorithm. Through 
this subroutine, the user may examine the stabilizability and detectability 
conditions implicit in the infinite-duration LQG methodology and, indirectly, 
compute minimal order state space realizations for transfer matrices. 

Finally, the analysis category of ORACLS includes subroutines (EXPSER, 
EXPADE, and EXPINT) for computing 
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and 


eAt 



dt 


These expressions are used in computing the transient response of linear time- 
invariant dynamical systems (TRANSIT). 

Subroutines are presented in the control law design part of ORACLS to 
implement some of the more common techniques of time-invariant LQG methodology. 
For the finite-duration optimal linear regulator problem with noise-free mea- 
surements, continuous dynamics, and integral performance index, a subroutine 
(CNTNREG) is provided to implement the negative exponential method for finding 
both the transient and steady-state solutions to the matrix Riccati equation 
(ref. 8). For the discrete version of this problem, the method of backward 
differencing is applied to find the transient and steady-state solutions to the 
discrete Riccati equation (DISCREG). For the infinite-duration optimal linear 
regulator problem with noise-free measurements, a subroutine is also included 
to solve the steady-state Riccati equation by the Newton algorithms described 
by Kleinman (ref. 9) for continuous problems and by Hewer (ref. 10) for dis- 
crete problems (RICTNWT). The methods described by Armstrong (ref. 11) and 
Armstrong and Rublein (ref. 12) are used to compute a stabilizing gain to ini- 
tialize the continuous and discrete Newton iterations (CSTAB and DSTAB) . A 
subroutine (PREFIL) is available for finding the prefilter gain to eliminate 
control-state cross-product terms in the quadratic performance index and another 
(SAMPL) computes the weighting matrices for the sampled-data optimal linear 
regulator problem (ref. 13). For cases with measurement noise, duality theory 
(ref. 4) and the foregoing optimal regulator algorithms are used to produce 
solutions to the continuous and discrete Kalman-Bucy filter problems (ASYMFIL). 
Finally, subroutines are included to implement the continuous (ref. 14) and 
discrete (ref. 15) forms of explicit (model in the system) and implicit (model 
in the performance index) model-following theory (EXPMDFL and IMPMDFL) . These 
subroutines generate linear control laws which cause the output of a time- 
invariant dynamical system to track the output of a prescribed model. 

In addition to the foregoing 43 primary purpose subroutines, ORACLS con- 
tains another 17 supporting subroutines used predominately by the algebraic 
equation and eigenvalue algorithms. All subroutines of the ORACLS package are 
(1) original codes by the author or (2) direct or author-modified copies of 
programs obtained from the VASP program (ref. 16), the software library of the 
Analysis and Computation Division at the Langley Research Center (LaRC), or the 
numerical analysis literature. The programs coded by the author and those 
obtained from the VASP program employ the dynamic storage capability in which 
data arrays are treated as packed one-dimensional arrays. Other subroutines 
were left in their original format. 

A detailed description of each ORACLS subroutine can be found in the 
following sections. In the last section, numerical examples are presented 
to illustrate the capability of ORACLS in both digital and continuous LQG con- 
troller design and, additionally, to demonstrate the construction of typical 
executive programs . 
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PRIMARY SUBROUTINES FOR INPUT-OUTPUT 


Subroutine RDTITL 

Description ; The purpose of RDTITL is to read a single card of Hollerith input 
which is loaded into the array TITLE of COMMON block LINES of RDTITL and auto- 
matically printed at the top of each page of output through the subroutine 
LNCNT. The Hollerith input is typically used to define, for future reference, 
the problem being solved by ORACLS. 

Subroutine RDTITL also serves to define certain data blocks important to other 
subroutines within ORACLS, such as information for COMMON blocks LINES and 
FORM discussed in the description of LNCNT and PRNT, respectively. The sub- 
routines BARSTW and SNVDEC require input parameters (EPSA, EPSB, and lAC) 
designating the accuracy to which solutions are to be obtained. When these 
programs are used internally to other subroutines, the accuracy parameters 
are set to values EPSAM, EPSBM, and lACM defined by DATA statements and con- 
tained in COMMON block TOL of RDTITL. Also, convergence parameters for ter- 
minating the recursive computations in subroutines SUM, EXPSER, EXPINT, SAMPL, 
DISCREG, CNTNREG, and RICTNWT are internally set. Parameters SUMCV (for SUM), 
RICTCV (for DISCREG, CNTNREG, and RICTNWT), MAXSUM (for SUM), and SERCV (for 
EXPSER, EXPINT, and SAMPL) are defined by DATA statements and contained in 
the COMMON block CONV. The user should specify the parameters of all COMMON 
blocks of RDTITL on the basis of his particular computing installation and 
problem to be solved. 

Subroutine RDTITL must be a part of every executive program provided by the 
user of ORACLS. If not, extraneous data appearing in the array TITLE will 
be printed at the top of each page of output and other COMMON block data will 
be undefined. 

Source of software : Ernest S. Armstrong, LaRC 


Calling sequence: 

CALL 

RDTITL 


Input arguments: 

None 



Output arguments: 

None 



COMMON blocks: LINES, 1 

FORM, TOL, 

CONV 

Error messages: 

None 



Field length: 32 

octal 

words (26 

decimal) 


Subroutine employed by RDTITL ; LNCNT 
Subroutines employing RDTITL; None 


Comments; None 
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Subroutine LNCNT 


Description: The purpose of LNCNT is to keep track of the number of lines 

printed and automatically paginate the output. Page length is controlled by 
the variable NLP set in the COMMON block LINES of subroutine RDTITL to 4^1. 
The variable NLP is an installation-dependent variable and may be changed as 
necessary. Subroutine LNCNT provides one line of print at the top of each 
page. This line contains 100 characters of which the first 80 are available 
for the programer's use and may be loaded by use of the subroutine RDTITL. 
The remainder contain "ORACLS PROGRAM.” The 100 characters are contained in 
the array TITLE within RDTITL. 

Source of software; VASP (ref. 16) with modifications by Ernest S. Armstrong, 
LaRC 

Calling sequence ; CALL LNCNT (N) 

Input argument ; 

N Number of lines to be printed 

Output arguments ; None 
COMMON block ; LINES 
Error messages ; None 

Field length ; 27 octal words (23 decimal) 

Subroutines employed by LNCNT ; None 

Subroutines employing LNCNT: RDTITL, PRNT, EQUATE, TRANP, SCALE, UNITY, TRCE, 

ADD, SUBT, MULT, JUXTC, JUXTR, FACTOR, SUM, BILIN, BARSTW, TESTSTA, EXPSER, 
EXPINT, VARANCE, CTROL, TRANSIT, SAMPL, PREFIL, CSTAB, DSTAB, DISCREG, 
CNTNREG, RICTNWT, ASYMREG, ASYMFIL, EXPMDFL, IMPMDFL, READ1 

Comments ; Subroutine LNCNT is completely internal to the ORACLS subroutines 
and the user need not refer to it unless he has a WRITE statement of his own. 
In that case, the user should put the statement CALL LNCNT(N) before each 
WRITE statement. 
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Subroutine READ 


Description ; The purpose of READ is to read from one to five matrices from 
cards along with their names and dimensions and print the same information. 
For each matrix, a header card is first read containing a four-character 
title followed by two integers giving the row and column size of the matrix 
using format (4H,4X,2l4). Then, the matrix data are read by rows using sub- 
routine READ1 (each row of the matrix starting on a new card) using format 
(8F10.2). Each matrix is then automatically printed using subroutine PRNT 
called from READ1 and packed by columns into one-dimensional arrays with the 
same names. 

Source of software; VASP (ref. 16) with modifications by Ernest S. Armstrong, 
LaRC 

Calling sequence; CALL READ(I,A,NA,B,NB,C,NC,D,ND,E,NE) 


Integer from 1 to 5 indicating the number of matrices to 
be read. For I < 5, the entries past I in the argument 
list may be omitted. 

Input matrices 


Two-dimensional vectors giving the number of rows and columns 
of the respective matrices and input through the header 
card; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 

If 0 is loaded for the row and column size, then the current 
matrix storage is unchanged, no data cards are read, and 
the previously stored matrix is printed. 

Output arguments : None 

COMMON blocks : None 

Error message : None directly from READ 

Field length : 155 octal words (109 decimal) 

Subroutine employed by READ ; READ1 
Subroutines employing READ; None 


Input arguments ; 
I 

A , B , C , 

D, E 

NA, NB, NC, 

ND, NE 


Comments: None 
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Subroutine PENT 


Description: The purpose of PENT is to print a single matrix with or without 

a descriptive heading either on the same page or on a new page. The descrip- 
tive heading, if desired, is printed before each matrix and is of the form 
" NAM MATEIX NA(1 ) EOWS NA(2) COLUMNS." The matrix is next printed by rows 
using format (1P7E16.7) for the first line and (3X1P7E16.7) for all subse- 
quent lines. Format for the printing is stored in COMMON block FOEM of 
EDTITL. 

Source of software: VASP (ref. 16) with modifications by Ernest S. Armstrong, 

LaEC 

Calling sequenc e: CALL PBNT(A,NA,NAM, lOP) 

Input arguments : 

A Matrix packed by columns in a one-dimensional array 

NA Two-dimensional vector giving the number of rows and columns of 

the matrix A: 

NA ( 1 ) = Number of rows of A 
NA(2) = Number of columns of A 

NAM Hollerith characters giving matrix name. Generally, NAM should 

contain 4 Hollerith characters and be written in the argument 
list as 4HXXXX. Alternatively, if 0 is inserted in the argument 
list for NAM, a blank name is printed. 

lOP Scalar print control parameter: 

1 Print heading and matrix on same page. 

2 Print heading and matrix after skipping to next page. 

3 Print only matrix with no heading on same page. 

4 Print only matrix with no heading after skipping to next 

page. 

Output arguments : None 

COMMON blocks : FOEM, LINES 

Error message: If NA(1) x NA(2) < 1 or NA(1) < 1, the message "EEBOB IN 

PENT MATEIX HAS NA = " is printed, and the program 

is returned to the calling point. 

Field length: 262 octal words (178 digital) 
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Subroutine employed by PRNT; LNCNT 


Subroutines employing PRNT: FACTOR, SUM, BILIN, BARSTW, TESTSTA, EXPSER, 

EXPINT, VARANCE, CTROL, TRANSIT, SAMPL, PREFIL, CSTAB, DSTAB, DISCREG, 
CNTNREG, RICTNWT, ASYMREG, ASYMFIL, EXPMDFL, IMPMDFL, READ1 


Comments : 


None 



PRIMARY SUBROUTINES FOR VECTOR-MATRIX OPERATIONS 


Subroutine EQUATE 

Description ; The purpose of EQUATE is to store a matrix in an alternate com- 
puter location. 

Source of software: VASP (ref. 16) with modifications by Ernest S. Armstrong, 

LaRC 

Calling sequence; CALL EQUATE(A,NA,B,NB) 


Input arguments ; 

A Matrix packed by columns in a one-dimensional array; not destroyed 

upon return 

NA Two-dimensional vector giving the number of rows and columns of A: 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

Output arguments ; 

B Matrix packed by columns in a one-dimensional array. Upon normal 

return, B = A. 

NB Two-dimensional vector: NB = NA upon normal return 

COMMON blocks ; None 

Error message ; If NA(1) x NA(2) <1 or NA(1 ) < 1, the message "DIMENSION 

ERROR IN EQUATE NA = •' is printed, and the program is returned 

to the calling point. 

Field length ; 46 octal words (38 decimal) 

Subroutine employed by EQUATE ; LNCNT 

Subroutines employing EQUATE: FACTOR, SUM, BILIN, BARSTW, TESTSTA, EXPSER, 

EXPINT, VARANCE, CTROL, TRANSIT, SAMPL, PREFIL, CSTAB, DSTAB, DISCREG, 
CNTNREG, RICTNWT, ASYMREG, ASYMFIL, EXPMDFL, IMPMDFL 

Comments: None 
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Subroutine TRANP 


Description ; The purpose of TRANP is to compute the transpose A* of a given 
matrix A. 

Source of software ; VASP (ref. 16) with modifications by Ernest S. Armstrong, 
LaRC 

Calling sequence ; CALL TRANP(A,NA,B,NB) 

Input arguments ; 

A Matrix packed by columns in a one-dimensional array; not destroyed 

upon return 

NA Two-dimensional vector giving the number of rows and columns of A; 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

Output arguments ; 

B Matrix packed by columns in a one-dimensional array. Upon normal 

return, B = A*. 

NB Two-dimensional vector: upon normal return, 

NB(1) = NA(2) 

NB(2) = NA(1) 

COMMON blocks ; None 

Error message : If NA(1) x NA(2) < 1 or NA(1) < 1, the message ’’DIMENSION 

ERROR IN TRANP NA = ” is printed, and the program is returned 

to the calling point. 

Field length ; 106 octal words (70 decimal) 

Subroutine employed by TRANP ; LNCNT 

Subroutines employing TRANP; FACTOR, SUM, BILIN, BARSTW, VARANCE, CTROL, 
TRANSIT, SAMPL, PREFIL, CSTAB, DSTAB, DISCREG, CNTNREG, RICTNWT, ASYMREG, 
ASYMFIL, EXPMDFL, IMPMDFL 

Comments; None 
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Subroutine SCALE 


Description ; The purpose of SCALE is to perform scalar multiplication on a 
given matrix. 

Source of software; VASP (ref. 16) with modifications by Ernest S. Armstrong, 
LaRC • 

Calling sequence ; CALL SCALE(A,NA,B,NB,S) 

Input arguments ; 

A Matrix packed by columns in one-dimensional array; not destroyed 

upon return 

NA Two-dimensional vector giving the number of rows and columns of A 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 

Not destroyed upon return 

S Scalar 

Output arguments ; 

B Matrix packed by columns in a one-dimensional array. Upon normal 

return, B = SA. 

NB Two-dimensional vector; NB = NA upon normal return 

COMMON blocks ; None 

Error message ; If NA(1) x NA(2) <1 or NA(1) < 1, the message "DIMENSION 

ERROR IN SCALE NA = " is printed, and the program is returned 

to the calling point. 

Field length ; 47 octal words (39 decimal) 

Subroutine employed by SCALE ; LNCNT 

Subroutines employing SCALE; FACTOR, BILIN, EXPSER, EXPINT, VARANCE, TRANSIT, 
SAMPL, PREFIL, CSTAB', DSTAB, CNTNREG, RICTNWT, ASYMREG, EXPMDFL, IMPMDFL 

Comments; None 
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Subroutine UNITY 


Desoription ; The purpose of UNITY is to generate an identity matrix. 

Source of software: VASP (ref. 16) with modifications by Ernest S. Armstrong, 

LaRC 

Calling sequence : CALL UNITY(A,NA) 

Input argument : 

NA Two-dimensional vector giving dimension of square identity matrix: 

NA(1) = NA(2) = Number of rows of A 
Not destroyed upon return 

Output argument : 

A Matrix packed by columns in a one-dimensional array. Upon normal 

return, A = I where I is an identity matrix of order NA(1 ) . 

COMMON blocks : None 

Error message : If NA(1) NA(2), the message "DIMENSION ERROR IN UNITY 

NA B " is printed, and the program is returned to the calling 

point. 

Field length : 6l octal words (49 decimal) 

Suly;;outine e mployed by UNITY : LMCNT 

Subroutine s employing UNITY : BILIN, EXPSER, EXPINT, TRANSIT 

Comments: None 
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Subroutine NULL 


Description ; The purpose of NULL is to generate a null matrix. 
Source of softwar e; Ernest S. Armstrong, LaRC 
Calling sequence ; CALL NULL(A,NA) 


Input argument ; 

NA Two-dimensional vector giving the momber of rows and columns of 

the desired null matrix; 

NA( 1 ) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

Output argument ; 

A Matrix packed by columns in a one-dimensional array. Upon 

normal return, A n 0 where 0 is a null matrix of order 
NA(1) X NA(2). 

COMMON blocks ; None 

Error message ; If NA(1) x NA(2) < 1 or NA(1) < 1, the message ’’DIMENSION 

ERROR IN NULL NA = " is printed, and the program is returned 

to the calling point. 

Field length ; 35 octal words (29 decimal) 

Subroutine employed by NULL ; LNCNT 

Subroutines employ ing NULL; BARSTW, CNTNREG 


Comments; 


None 



Subroutine TRCE 


Desoription ; The purpose of TRCE is to compute the trace of a square matrix. 

Source of software ; VASP (ref. 16) with modifications by Ernest S. Armstrong, 
LaRC 

Calling sequence ; CALL TRCE( A,NA,TR) 

Input arguments ; 

A Matrix packed by columns in a one-dimensional array; not destroyed 

upon return 

NA Two-dimensional vector giving number of rows and columns of A; 

NA(1) = NA(2) = Order of A 
Not destroyed upon return 

Output argument ; 

TR Scalar. Upon normal return, TR = Trace of A. 

COMMON blocks ; None 

Error message ; If NA(1) NA(2), the message "TRACE REQUIRES SQUARE MATRIX 

NA = " is printed, and the program is returned to the calling 

point . 

Field length ; 52 octal words (42 decimal) 

Subroutine employed by TRCE ; LNCNT 
Subroutine employing TRCE ; EXPSER 
Comments : None 


\ 
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Subroutine ADD 


Desoription ; The purpose of ADD is to perform matrix addition C = A + B for 
given matrices A and B. 

Source of software: VASP (ref. 16) with modifications by Ernest S. Armstrong, 

LaRC 

Calling sequence ; CALL ADD(A,NA,B,NB,C,NC) 

Input arguments ; 

A, B Matrices packed by columns in one-dimensional arrays; not destroyed 

upon return 

NA, NB Two-dimensional vectors giving number of rows and columns of respec- 
tive matrices; for example, 

NA( 1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

Output arguments ; 

C Matrix packed by columns in a one-dimensional array. Upon normal 

return, C = A + B. 

NC Two-dimensional vector: upon normal return, 

NC(1) = NA(1) 

NC(2) = NA(2) 

COMMON blocks : None 

Error message: If either NA(1) ^ NB(1), NA(2) ^ NB(2), NA(1 ) < 1 , or 

NA(1) X NA(2) < 1, the message "DIMENSION ERROR IN ADD NA = , 

NB = " is printed, and the program is returned to the calling 

point . 

Field length : 56 octal words (46 decimal) 

Subroutine employed by ADD : LNCNT 

Subroutines employing ADD: SUM, EXPSER, EXPINT, TRANSIT, SAMPL, DSTAB, DISCREG, 

CNTNREG, RICTNWT, ASYMREG, EXPMDFL, IMPMDFL 

Comments: None 
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Subroutine SUBT 


Description ; The purpose of SUBT is to perform matrix subtraction C = A - B 
for given matrices A and B. 

Source of software: VASP (ref. 16) with modifications by Ernest S. Armstrong, 

LaRC 

Calling sequence: CALL SUBT(A,NA,B,NB,C,NC) 


Input arguments ; 

A, B Matrices packed by columns in one-dimensional arrays; not destroyed 

upon return 

NA, NB Two-dimensional vectors giving the number of rows and columns of 
respective matrices; for example, 

NA( 1 ) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

Output arguments ; 

C Matrix packed by columns in a one-dimensional array. Upon normal 

return, C = A - B. 

NC Two-dimensional vector: upon normal return, 

NC(1) = NA(1) 

NC(2) = NA(2) 

COMMON blocks ; None 

Error message: If either NA(1) NB( 1 ) , NA(2) ^ NB(2), NA(1) < 1 , or 

NA(1) X NA(2) < 1, the message "DIMENSION ERROR IN SUBT NA = , 

NB = " is printed, and the program is returned to the calling 

point. 

Field length ; 56 octal words (46 digital) 

Subroutine employed by SUBT ; LNCNT 

Subroutines employing SUBT: TRANSIT, PREFIL, CSTAB, DSTAB, DISCREG, CNTNREG, 

RICTNWT, ASYMREG, EXPMDFL, IMPMDFL 

Comments: None 
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Subroutine MULT 


Description ; The purpose of MULT is to perform matrix multiplication C s AB 
for given matrices A and B. 

Source of software ; VASP (ref. 16) with modifications by Ernest S. Armstrong, 
LaRC 

Calling sequence; CALL MULT(A,NA,B,NB,C,NC) 


Input arguments ; 

A, B Matrices packed by columns in one-dimensional arrays; not destroyed 

upon return 

NA, NB Two-dimensional vectors giving the number of rows and columns of 
respective matrices; for example, 

NA( 1 ) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

Output arguments ; 

C Matrix packed by columns in a one-dimensional array. Upon normal 

return, C = AB. 

NC Two-dimensional vector; upon normal return, 

NC(1) = NA(1) 

NC(2) = NB(2) 

COMMON blocks ; None 

Error message: If either NA(2) NB(1), NA(1) < 1, NA( 1) x NA(2) < 1 , or 

NA(1) X NB(2) < 1, the message "DIMENSION ERROR IN MULT NA = , 

NB = " is printed, and the program is returned to the calling 

point . 

Field length ; 154 octal words (108 decimal) 

Subroutine employed by MULT ; LNCNT 

Subroutines employing MULT: FACTOR, SUM, BILIN, EXPSER, EXPINT, VARANCE, CTROL, 

TRANSIT, SAMPL, PREFIL, CSTAB, DSTAB, DISCREG, CNTNREG, RICTNWT, ASYMREG, 
EXPMDFL, IMPMDFL 

Comments : None 
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Subroutine MAXEL 


Description : The purpose of MAXEL is to compute the maximum of the absolute 

values of the elements of a real matrix. 

Source of software : Ernest S. Armstrong, LaRC 

Calling sequence : CALL MAXEL( A,NA,ELMAX) 

Input arguments : 

A Matrix packed by columns in a one-dimensional array; not destroyed 

upon return 

NA Two-dimensional vector giving the number of rows and columns of A: 

NA( 1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

Output argument ; 

ELMAX Scalar. Upon normal return, ELMAX is the maximum of the absolute 
values of the elements of A. 

COMMON blocks : None 

Error messages : None 

Field length : 27 octal words (23 decimal) 

Sub routines employed by MAXEL ; None 

Subroutines employing MAXEL: SUM, EXPSER, EXPINT, SAMPL, DISCREG, CNTNREG, 

~ rictnW 

Comments: None 
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Subroutine NORMS 


Desoription ; The purpose of NORMS is to compute either the 5<i , JI2 

(Euclidean), or 5.^^ matrix norms for a real m x n matrix A stored as a 
variable-dimensioned two-dimensional array. The norms 9 , 1 , 9,2 > ^<x> 

defined, respectively, as 


m 

IIa|Ii = oAx } 

ISkSn 


hjk| 


IUI2 = 





1/2 



max 
1< j<m 


I 

U-1 


l^jkl 


Source of software; LaRC Analysis and Computation Division subprogram library 


Calling 

sequence: 

CALL NORMS ( MAXROW , M , N , A , lOPT , RLNORM) 

Input arguments: 


MAXROW 

Maximum first dimension of array A as given in the DIMENSION 


statement of the calling program 

M 

Number 

of rows of matrix A 

N 

Number 

of columns of the matrix A 

A 

Matrix 

whose norm is desired stored in a two-dimensional array 

lOPT 

Scalar 

norm selector: 


1 

Compute 9,-j. 


2 

Compute 9,2 . 


3 

Compute 9,g^. 

Output argument: 


RLNORM 

Scalar. 

Upon normal return, RLNORM is the appropriate norm. 

COMMON blocks : None 
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Error messages; None 


Field length ! 152 octal words (106 decimal) 

Su broutines employed by NORMS ; None 

Subroutines employing NORMS ; BILIN, EXPSER, EXPINT, SAMPL, CSTAB 

Comments : NORMS can be applied to matrices stored in packed one-dimensional 

arrays by placing MAXROW = M in the calling sequence. 
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Subroutine JUXTC 


Description ; The purpose of JUXTC is to construct a matrix [A,B] from given 
matrices A and B. 

Source of software: VASP (ref. 16) with modifications by Ernest S. Armstrong, 

LaRC 

Calling sequence; CALL JUXTC(A,NA,B,NB,C,NC) 


Input arguments ; 

A, B Matrices packed by columns in one-dimensional arrays; not destroyed 

upon return 

NA, MB Two-dimensional vectors giving the number of rows and columns of 
the respective matrices; for example, 

NA( 1 ) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

Output arguments : 

C Matrix packed by columns in a one-dimensional array. Upon normal 

return, C = [A,B]. 

NC Two-dimensional vector: upon normal return, 

NC(1) = NA(1) 

NC(2) = NA(2) + NB(2) 

COMMON blocks ; None 

Error message: If either NA( 1 ) ^ NB(1), NA( 1 ) < 1, NA(1) x NA(2) < 1 , or 

NA(2) + NB(2) < 1, the message "DIMENSION ERROR IN JUXTC NA = , 

NB = " is printed, and the program is returned to the calling 

point . 

Field length ; 76 octal words (62 digital) 

Subroutine employed by JUXTC ; LNCNT 

Subroutines employing JUXTC; TESTSTA, CTROL, DSTAB, ASYMREG 


Comments : 


None 


Subroutine JUXTR 


Description ; The purpose 

matrices A and B. 

Source of software: VASP 

LaRC 


of JUXTR is to construct a matrix 



from given 


(ref. 16) with modifications by Ernest S. Armstrong, 


Calling sequence ; CALL JUXTR(A,NA,B,NB,C,NC) 

Input arguments ; 

A, B Matrices packed by columns in one- dimensional arrays; not destroyed 

upon return 

NA, NB Two-dimensional vectors giving the number of rows and columns of 
the respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 


Output arguments; 


Matrix packed by columns in a one- dimensional array. 

‘a‘ 


Upon normal 


return, C = 


B 


NC Two-dimensional vector: upon normal return, 

NC(1) = NA(1) + NB(1) 

NC(2) = NA(2) 


COMMON blocks : None 

Error message: If either NA(2) NB(2), NA(1) < 1, NA( 1 ) x NA(2) < 1 , or 

NA(2) < 1, the message "DIMENSION ERROR IN JUXTR NA = , 

NB = ” is printed, and the program is returned to the calling 

point . 

Field length ; 144 octal words (100 digital) 

Subroutine employed by JUXTR ; LNCNT 
Subroutines employing JUXTR; BARSTW, CNTNREG 


Comments : None 
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PRIMARY SUBROUTINES FOR ANALYSIS OF CONSTANT LINEAR SYSTEMS 


Subroutine FACTOR 

^Description : The purpose of FACTOR is to compute a real m x n (m ^ n) 

matrix D of rank m such that a real n x n nonnegative definite 
matrix Q can be factored as 

Q = D'D 

The method is first to perform a singular- value (eigenvalue since Q = Q* SO) 
decomposition of Q as 

Q = PJP' 


where J is a diagonal matrix of eigenvalues of Q, and then define D as 


{j P' 


after eliminating those rows of '/T with all zero elements. Eigenvalues 
of Q are considered negligible (zero) if they are less than ||q|| x tO“(IAC) 
using the matrix Euclidean norm. 

Source of software : Ernest S. Armstrong, LaRC 

Calling sequence: CALL FACTOR ( Q, NQ, D, ND, lOP ,IAC, DUMMY) 


Input arguments : 

Q Nonnegative definite matrix packed by columns in a one-dimensional 

array; not destroyed upon return 

NQ Two-dimensional vector giving the number of rows and columns of Q: 

NQ(1) = NQ(2) = Number of rows of Q 
Not destroyed upon return 

lOP Scalar print parameter: 

0 Return within printing results. 

Otherwise Print Q, D, and D'D. 

lAC Scalar parameter to be used for zero test of eigenvalues of Q 

DUMMY Vector of working space for computations with dimension at least 

n^ + n 
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Output arguments; 


D Matrix packed by columns in a one-dimensional array with dimension 

at least n^. Upon normal return, D'D = Q within numerical 
accuracy. 

ND Two-dimensional vector giving number of rows and columns of D: 

upon normal return, 

ND(1) = m 
ND(2) = n 

COMMON blocks ; None 

Error messages ; > 

(1) If the singular-value decomposition subroutine SNVDEC fails to converge 

to a singular value after 30 iterations, the message "IN FACTOR, SNVDEC 

HAS FAILED TO CONVERGE TO THE SINGULAR VALUE AFTER 30 ITERATIONS" 

is printed, and the program is returned to the calling point, 

(2) If an eigenvalue of Q is greater than but close to the value 

ZTEST = llQll X lO-(IAC) (less than 16 x ZTEST), the message "IN FACTOR, 
THE MATRIX Q SUBMITTED TO SNVDEC IS CLOSE TO A MATRIX OF LOWER RANK 

USING ZTEST a / IF THE ACCURACY IS REDUCED THE RANK MAY ALSO BE 

REDUCED / CURRENT RANK = " is printed along with the computed 

singular values, and the computation continues. 

Field length ; 422 octal words (274 decimal) 

Subroutines employed by FACTOR ; EQUATE, SNVDEC, LNCNT, TRANP, PRNT, MULT 

Subroutines employing FACTOR; None 


Comments: None 


25 


Subroutine EIGEN 


Description ; The purpose of EIGEN is to compute all the eigenvalues and 
selected eigenvectors of a real n x n matrix A stored as a variable- 
dimensioned two-dimensional array. The input matrix is first balanced by 
exact similarity transformations such that the norms of corresponding rows 
and columns are nearly equal (ref. 17). The balanced matrix is reduced to 
upper Hessenberg form by stabilized elementary similarity transformations 
(ref. 18). All of the eigenvalues of the Hessenberg matrix are found by the 
double shift QR algorithm (ref. 19). The desired eigenvectors of the 
Hessenberg matrix are then found by the inverse iteration method (ref. 2). 

Source of software ; LaRC Analysis and Computation Division subprogram library 
with modifications by Ernest S. Armstrong, LaRC 

T 

Calling sequence; CALL EIGEN(MAX,N,A,ER,EI,ISV,ILV,V,WK,IERR) 


Input arguments ; 

MAX Maximum first dimension of the array A as given in the DIMENSION 

statement of the calling program 

N Actual order of matrix (n) 

A Matrix whose eigenvalues and selected eigenvectors are desired 

stored in a real two-dimensional array. The contents of this 
array are destroyed upon return. 

ISV The number of eigenvalues, smallest in absolute value, for which 

eigenvectors are desired counting complex conjugates 

ILV The number of eigenvalues, largest in absolute value, for which 

eigenvectors are desired counting complex conjugates 

WK Vector of working space for computations of dimension; 

3N if ISV + ILV = 0 

N(N+7) otherwise 

Output arguments ; 

ER One-dimensional real array containing the real parts of the eigen- 

values, dimensioned at least N in the calling program 

El One -dimensional real array containing the imaginary parts of the 

eigenvalues, dimensioned at least N in the calling program 

ISV On output, ISV is the number of eigenvalues, smallest in absolute 

value, for which eigenvectors were computed counting complex 
conjugates 
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ILV On output, ILV is the number of eigenvalues, largest in absolute 

value, for which eigenvectors were computed counting complex 
conjugates 

V Two-dimensional array containing the eigenvectors, normalized to 

unit length on a normal return. It suffices to have V dimen- 
sioned MAX as first dimension and N as second. 

lERR Integer error code: 

lERR = 0 Normal return 

lERR = -J Vectors for the Jth eigenvalue did not converge to 
an eigenvector. Appropriate column (columns if 
Jth eigenvalue is complex) of V is set to zero. 

If failure occurs more than once, the index for 
the last such occurrence is in lERR. 
lERR = J The Jth eigenvalue has not been determined after 
30 iterations of the QR algorithm. 

COMMON blocks : None 

Error messages : None. User should test lERR after return. 

Field length : 1131 octal words (601 decimal) 

Subroutines employed by EIGEN : BALANC, ELMHES, HQR, INVIT, ELMBAK, BALBAK 

Subroutines employing EIGEN : BILIN, TESTSTA, CSTAB, DSTAB, CNTNREG, ASYMREG 

Comments : Upon normal return, eigenvalues are stored in ascending magnitude 

with complex conjugates stored with positive imaginary parts first. The 
eigenvectors are packed and stored in V in the same order as their eigen- 
values appear in ER and El. Only one eigenvector is computed for complex 
conjugates (for conjugate with positive imaginary part). Upon error exit -J, 
eigenvalues are correct and eigenvectors are correct for all nonzero vectors. 
Upon error exit J, eigenvalues are correct but unordered for indices 
IERR+1 ,IERR+2, . . . ,N, and no eigenvectors are computed. 

EIGEN may be used for matrices packed as a one-dimensional array by setting 
MAX = N in the calling sequence. If all eigenvectors are desired, set 
ISV = N and ILV = 0. 
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Subroutine SYMPDS 


Description ; The purpose of SYMPDS is to solve the matrix equation, 

AX = B 

where A is a symmetric positive definite matrix and B is a matrix of 
constant vectors. The determinant of A may also be evaluated. Solution 
is by Cholesky decomposition (ref. 2): A is factored as As LDL', where 

L is a unit lower triangular matrix and D is a positive definite diagonal 
matrix. An option is provided for computing only the Cholesky factorization 
of A without solving a complete matrix equation. Both A and B are 
stored as variable-dimensioned two-dimensional arrays. 

Source of software ; LaRC Analysis and Computation Division subprogram library 
with modifications by Ernest S. Armstrong, LaRC 

Calling sequence; CALL SYMPDS (MAXN,N, A, NRHS,B,IOPT, IF AC, DETERM, ISCALE,P,IERR) 


Input arguments ; 

MAXN The maximum first dimension of A as given in the DIMENSION 

statement of the calling program 

N The number of rows of A 

A Coefficient matrix stored as a variable-dimensioned two-dimensional 

array. Two types of input are possible; the first is the unde- 
composed coefficient matrix, and the second is the Cholesky 
decomposition A = LDL*. For A in unfactored form, the con- 
tents are destroyed upon return. If the factored form of A is 
input, A contains the upper triangular elements of the coeffi- 
cient matrix and the elements of the unit lower triangular matrix 
L except the diagonal elements of L which are understood to be 
all unity. The reciprocals of the elements of D are input in 
the array P . 

NRHS The number of column vectors of the matrix B. No data are required 

if only factorization of A is desired, but NRHS still must 
appear as an argument of the calling sequencfe. 

B Two-dimensional array that must have first dimension MAXN and 

second dimension at least NRHS in the calling program. On 
input, B contains the elements of the constant vectors, and B 
is destroyed upon return. If only a factorization of A is 
required, no data are necessary, but B still must appear as an 
argument of the calling sequence. 

lOPT Integer for determinant evaluation; 

0 Determinant is not evaluated. 

1 Determinant is evaluated. 
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IFAC Integer specifying whether or not the Cholesky decomposition of the 

coefficient matrix is to be computed: 

0 Cholesky decomposition for the matrix A is to be computed ; 

and equation solved. 

1 Cholesky decomposition form of A is input so no decomposi- 

tion is required. 

2 Only Cholesky decomposition of A is required. 

P One-dimensional array dimensioned at least N in the calling pro- 

gram. If the unfactored form of A is input, nothing need be 
input for P. If the factored form of A is input, P contains 
the reciprocals of the diagonal elements of D. 

Output arguments : 

A Upon normal return, the array A contains the original elements 

of the matrix A and the elements of the unit lower triangular 
matrix L except the diagonal elements of L which are under- 
stood to be all unity. 

B Upon normal return, if a system of equations is to be solved, each 

solution vector of X is stored over the corresponding constant 
vector of the input array B. 

DETERM, Determinant evaluation parameters; 

ISCALE det (A) = DETERM x iq( lOOxiSCALE) 

P Upon normal return, P contains the reciprocals of the diagonal 

elements of D. 

lERR Error code : 

0 Normal return 

1 A is not symmetric positive definite. 

COMMON blocks ; None 

Error messages ; None. The parameter lERR should be tested after return. 

Field length ; 373 octal words (251 decimal) 

Subroutines employed by SYMPDS ; None 

Subroutines employing SYMPDS ; PREFIL, CNTNREG, RICTNWT, EXPMDFL 

Comments ; SYMPDS may be used for matrices packed as a one-dimensional array 
by setting MAXN = N in the calling sequence. 
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Subroutine GELIM 


Description; The purpose of GELIM is to solve the real matrix equation, 


AX = B 


where A is required to be square and nonsingular and B is a matrix of 
constant vectors. Solution is by Gaussian elimination or LU factorization 
(ref. 2) in which A is factored as 


PA = LU 


where L is a unit lower triangular matrix, U is an upper triangular 
matrix, and P is a permutation matrix representing the row pivotal strategy 
associated with the LU factorization. Both A and B are stored as 
variable-dimensioned two-dimensional arrays. 

Source of software ; LaRC Analysis and Computation Division subprogram library 

Calling sequence; CALL GELIM(NMAX,N,A,NRHS,B,IPIVOT,IFAC,WK,IERR) 


Input arguments : 

NMAX The maximum first dimension of the array A as given in the 

DIMENSION statement of the calling program 

N The number of rows of A 

A Coefficient matrix stored as a variable-dimensioned two-dimensional 

array. Two types of input are possible: the first is the unfac- 

tored coefficient matrix, and the second is the triangular fac- 
torization A = LU. For A input in factored form, A = (L\U) 
should be used neglecting the unity elements of L, and the 
pivotal strategy employed should be input through the array 
IPIVOT. For A input in unfactored form, input data are 
destroyed . 

NRHS Number of column vectors of the matrix B 

B Two-dimensional array that must have first dimension NMAX and 

second dimension at least NRHS in the calling program. On input, 
B contains the elements of the constant vectors and is destroyed 
upon return. 
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IPIVOT An integer array dimensioned at least N in the calling program. 

If the factored form of A is input, IPIVOT contains the pivotal 
strategy by the rule, 


IPIVOT(I) = J 


which states that row J of matrix A was used to pivot for the 
Ith unknown. 

IFAC Factorization parameter: 

0 Compute L, U, and pivotal strategy. 

1 Do not compute L, U, and pivotal strategy; the factoriza- 

tion and strategy are input. 

WK A one-dimensional array dimensioned at least N in the calling 

program and used as a work storage array 

Output arguments : 

A Upon normal return, the unit lower and upper matrices are over- 

stored in A as A = (L\U) neglecting the unity elements L. 

B Upon normal return, each solution vector of X is stored over 

the corresponding constant vector of the input array B. 

IPIVOT Upon normal return, IPIVOT contains the pivotal strategy as 
previously explained. 

lERR Singularity test parameter: 

0 A is nonsingular. 

1 A is singular. 

COMMON blocks ; None 

Error messages ; None. Upon return the parameter lERR should be tested. 

Field length ; 261 octal words (177 digital) 

Subroutine employed by GELIM ; DETFAC 

Subroutines employing GELIM ; BILIN, TRANSIT, DSTAB, CNTNREG 

Comments : If an LU factorization of A is desired without a complete equa- 
tion solution, the subroutine DETFAC may be employed. 

GELIM may be used for matrices packed as a one-dimensional array by setting 
NMAX = N in the calling sequence. 
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Subroutine SNVDEC 


Desoription ; The purpose of SNVDEC is to compute the singular-value decomposi- 
tion (ref. 2) of a real m x n (m S n) matrix A by performing the factori- 
zation, 


A = UQV 


where U is an m x n matrix whose columns are n orthonormalized eigen- 
vectors associated with the n largest eigenvalues of AA', V is an 
n X n matrix whose columns are the orthonormalized eigenvectors associated 
with the n eigenvalues of A 'A, and 

Q = diag (ai , 02 , • •• 

where (i □ 1,2,...,n) are the nonnegative square roots of the eigen- 
values of A' A, called the singular values of A. Options are provided for 
the computation of rank A, singular values of A, an orthonormal basis for 
the null space of A, the pseudoinverse of A, and the least squares solution 
to 


AX = B 


Both A and B are stored as variable-dimensioned two-dimensional arrays. 

The computational procedure is described in reference 2 on pages 135-151. 
Basically, Householder transformations are applied to reduce A to bidiagonal 
form after which a QR algorithm is used to find the singular values of the 
reduced matrix. Combining results gives the required construction. 

Source of software ; LaRC Analysis and Computation Division subprogram library 
with modifications by Ernest S. Armstrong, LaRC 

Calling sequence : CALL SNVDEC ( lOP , MD , ND , M , N , A , NOS , B , lAC , ZTEST , Q , V , IRANK , 

APLUS,IERR) 

Input arguments ; 

lOP Option code : 

1 The rank and singular values of A will be returned. 

2 The matrices U and V will be returned in addition to the 

information for lOP = 1 . 

3 In addition to the information for lOP = 2, the least 

squares solution to AX = B will be returned. 

H The pseudoinverse of A will be returned in addition to the 
information for lOP = 2. 

5 The least squares solution will be returned in addition to 
the information for lOP = 4. 
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MD The maximum first dimension of the array A as given in the 

DIMENSION statement of the calling program 

ND Maximum first dimension of the array V 

M The number of rows of A 

N The number of columns of A 

A Matrix stored as a variable-dimensioned two-dimensional array. 

Input A is destroyed. 

NOS The number of column vectors of the matrix B 

B Two-dimensional array that must have row dimension at least NOS 

in the calling program. B contains the right sides of the 
equation to be solved for lOP =3 or lOP =5* B need not 
be input for other options but must appear in the calling 
sequence . 

lAC The number of decimal digits of accuracy in the elements of the 

matrix A. This parameter is used in the test to determine zero 
singular values and thereby the rank of A. 

Output arguments ; 

A On normal return, A contains the orthogonal matrix U except 

when lOP = 1 . 

B On normal return, B contains the least squares solution for 

lOP =3 or lOP = 5. 

ZTEST The zero test computed as ||Al| x using the matrix 

Euclidean norm except when N = 1 . When N = 1 , 

ZTEST = 


Q A one-dimensional array of dimension at least N which upon return 

contains the singular values in descending order 

V A two-dimensional array that must have first dimension ND and 

second dimension at least N. Upon normal return, this array 
contains the orthogonal matrix V except when lOP □ 1 . The 
last N - IRANK columns of V form a basis for the null space 
of A. 

IRANK Rank of the matrix A determined as the number of nonzero singular 
values using ZTEST 
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APLUS 


lERR 


A two-dimensional array of first dimension ND and second dimension 
at least M. Upon normal return, this array contains the pseudo- 
inverse of the matrix A. If lOP does not equal 4 or 5, this 
array need not be dimensioned, but a dummy parameter must appear 
in the calling sequence. 

Error indicator; 

lERR =0 A normal return 

lERR = K > 0 The Kth singular value has not been found after 

30 iterations of the QR algorithm procedure. 
lERR = -1 Using the given lAC, A is close to a matrix 

which is of lower rank and if the accuracy is 
reduced, the rank of the matrix may also be 
reduced . 


COMMON blocks : None 

Error messages ; None. The user should examine lERR after return. 

Field length ; 2072 octal words (1082 decimal) 

Subroutines employed by SNVDEC ; None 

Subroutines employing SNVDEC ; FACTOR, CTROL, CSTAB, DSTAB, DISCREG 

Comments ; SNVDEC may be applied to matrices stored as one-dimensional arrays 
by setting MD = M and ND = N in the calling sequence. 

The subroutine is internally restrict.ed to N S 150. 


Subroutine SUM 


Description ; The purpose of SUM is to evaluate until convergence the matrix 
series 


00 



i=0 


where A and C are n x n and m x m real constant matrices, respec- 
tively. The matrix B is real constant and n x m. The series is numeri- 
cally summed by successively evaluating the partial sum sequence 


S(j+1) = S(j) + U(j) S(j) V(j) 
U(j+1) = U2(j) > 

V(j+1) = V2(j) ^ 

with 

S(0) = B 


(j = 0,1,...) 


U(0) = A 


V(0) = C 


The symbol S(j) represents the sum of the first 2J terms of the power 
series. Evaluation of the sequence continues until the number of terms evalu- 
ated exceeds MAXSUM, which is specified in the COMMON block CONV of subroutine 
RDTITL, or until convergence is reached. Numerically, convergence of the 
S(j) sequence is determined by testing the improvement in the element of 
S(j) of largest magnitude (measured relatively if the magnitude is less than 
unity, and absolutely otherwise) . The sequence is assumed to have converged 
when this improvement is less than the parameter SUMCV also found in the 
COMMON block CONV of subroutine RDTITL. A condition sufficient for the con- 
vergence of the power series, independent of B, is that each eigenvalue 
of A and C have complex modulus less than unity. When the series con- 
verges, it represents a solution X to 

X = AXC + B 


Source of software; Ernest S. Armstrong, LaRC 
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Calling sequence; CALL SUM(A,NA,B,NB,C,NC,IOP,SYM, DUMMY) 


Input arguments; 


A, B, C Compatible matrices packed by columns in one-dimensional 

arrays. A, B, and C are destroyed upon return. 

NA, NB, NC Two-dimensional vectors giving the number of rows and columns 
of the respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 


lOP Print parameter: 

0 Do not print results. 

Otherwise Print input A, B, C, and the sum denoted 
as X. 


SYM Logical variable: 

TRUE Indicates A = C 

FALSE Otherwise 

DUMMY Vector of working space for computations with dimension at 

least maximum of (n^ ,m^ ,2nm) 


Output argument: 


B Upon normal return, the sum X is stored in B. 

COMMON block : CONV 

Error message: If the number of terms in the partial sum sequence S(j) 

exceeds MAXSUM, the message "IN SUM, THE SEQUENCE OF PARTIAL SUMS HAS 
EXCEEDED STAGE WITHOUT CONVERGENCE" is printed. 

Field length : 367 octal words (2M7 decimal) 

Subroutines employed by SUM : PRNT, MULT, MAXEL, ADD, EQUATE, TRANP, LNCNT 

Subroutines employing SUM : BILIN, VARANCE, RICTNWT, EXPMDFL 

Comments : A method for scaling the A and C matrices in order to possibly 

improve convergence is found in reference 15. Also found in reference 15 is 
a bilinear transformation method for converting the equation. 


X = AXC + B 


into one solvable (assuming a unique solution exists) by the subroutine 
BARSTW of ORACLS. 
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Subroutine BILIN 


Description; The purpose of BILIN is to solve the matrix equation, 


AX + XB = C 


where A and B are real constant matrices of dimension n x n and m x m, 
respectively. The matrix C is real constant and of dimension n x m. It 
is assumed that all eigenvalues of A and B have strictly negative real 
parts. The method of solution employs the bilinear transformation technique 
described by Smith (ref. 6), wherein it is established that the foregoing 
X equation is equivalent to 


X = UXV + W 


with 


U = (3In - A)-‘'(6In + A) 

V = (3Im + B)(3Im - B)-1 
W = -23(3In - A)-1c(3Im - B)"’' 


where 3 is a scalar greater than zero and I^ and I^ are n x n and 
m X m identity matrices. The eigenvalues of U and V lie within the 
unit circle in the complex plane; therefore, the series 


I 

i -r 


Uiwvi 


converges and represents X. The subroutine SUM is used to evaluate the 
infinite series. 

The user has the option of inputing 3 externally or having the program com- 
pute 3 internally based on the eigenvalues of A and B as follows. Let 


Xj = aj + ibj ( j = 1 ,2, . . . ,n) 
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be eigenvalues of A ordered so that 


|Xil ^ IX 2 I ^ S |Xn| 

Choosing fi'om the equation, 

(3q + a-|)^ + b-]^ (Bq + ajj)^ + bjj^ 

(Bq - ai)2 + b-]2 (3 q _ ajj)2 + bjj2 

gives 

ai(an^ + bn^) - an(ai2 + bi^) 

Bo 2 = p ^ 

(an - ai) 

When a^ - a-| = 0 or Bq^ = 0> instead set 

^ (aj2 + bj2)^/2 

Bo 

n 

If the eigenvalue computation fails, set 
Bo = 2|1a1| 

using the 2-i matrix norm (see subroutine NORMS). Similarly, compute B-| 
based on B and put 

1 

B = -(Bo + Bi) 

Source of software : Ernest S. Armstrong, LaRC 

Calling sequence; CALL BILIN(A,NA,B,NB,C,NC, lOP , BETA, SYM, DUMMY) 
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Input arguments; 


A, B, C Compatible matrices packed by columns in one-dimensional arrays. 

C is destroyed upon return, but A and B are not. If 
A = B ' , no data need be entered for B and NB , but their 
symbols must still appear in the argument list. 

NA, NB, NC Two-dimensional vectors giving the number of rows and columns 
of the respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 


lOP 


BETA 


Two-dimensional parameter vector; 

I0P(1) □ 0 Do not print results. 

Otherwise Print A, B, C, BETA, and X. 


I0P(2) = 0 Use input value BETA for g. 
Otherwise Compute g as previously described. 


Scalar value g for numerical conditioning. No input is 
required if I0P(2) is nonzero. 


SYM Logical variable: 

TRUE Indicates A = B‘ 

FALSE Otherwise 


DUMMY 


Vector of working space for computations with dimension at 
least (4p^ + 2p) with p = max(n,m) 


Output argument; 


C Upon normal return, the solution X is stored in C. 

COMMON blocks ; None 
Error messages ; 

( 1 ) If g is computed internally and the eigenvalue computation for A 

fails, the message "IN BILIN, THE EIGENVALUE OF A HAS NOT BEEN 

DETERMINED AFTER 30 ITERATIONS" is printed. A similar message is 
printed if the eigenvalue computation for B fails. 

(2) If the (gin - A)"’’ computation fails because of singularity, the 

message "IN BILIN, THE MATRIX (BETA)I - A IS SINGULAR, INCREASE BETA" 
is printed. A similar message is printed if the (3Im “ B)“^ compu- 
tation fails. 
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Field length; 1512 octal words (.8H2 decimal) 


Subroutines employed by BILIN: LNCNT, PRNT, TRANP, EQUATE, EIGEN, SCALE, MULT, 

SUM, GELIM, NORMS 

Subroutines employing BILI N; VARANCE, CSTAB, RICTNWT 

Comments ; For A and B not satisfying the requirement of eigenvalues with 
strictly negative real parts, but still admitting a unique solution to 


AX + XB = C 


the subroutine BARSTW may be applied. 
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Subroutine BARSTW 


Description; The purpose of BARSTW is to solve the matrix equation 


AX + XB = C 


where A and B are real constant matrices of dimension n x n and m x m, 
respectively. The matrix C is real constant and of dimension n x m. It 
is assumed that 

A B 


A B 

where Xi and Xj are eigenvalues of A and B, respectively. The matrix 
equation then has a unique solution X. The method of solution is based on 
transforming A and B to real Schur form as described by Bartels and 
Stewart (ref. 5). 

Source of software ; Ernest S. Armstrong, LaRC 

Calling sequence ; CALL BARSTW(A,NA,B,NB,C,NC,IOP,SYM,EPSA,EPSB, DUMMY) 


Input arguments; 


A, B, C Compatible matrices packed by columns in one-dimensional arrays. 

C is destroyed upon return, but A and B are not. If 
A B B ' and C = C ' , no data need be entered for B and NB , 
but their symbols should appear in the argument list. 


NA, NB, NC Two-dimensional vectors giving the number of rows and columns 
of the respective matrices; for example, 

NA ( 1 ) B Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 


lOP 


SYM 


Print option parameter: 

0 Do not print. 

Otherwise Print A, B, C, and X. 

Logical variable: 

TRUE Indicates A = B' and C = C’ 

FALSE Otherwise 


EPSA, EPSB Convergence criteria for the reduction of A and B to Schur 
form. EPSA should be set slightly smaller than 10^"^^, 
where N is the number of significant digits in the elements 
of A. EPSB is similarly defined from B. 



DUMMY Vector of working space for computations with dimensions at 

least 

2(n + 1)2 for SYM = TRUE 

2(n + l-)2 + 2(m + 1)2 for SYM = FALSE 

Output argument ; 

C Upon normal return, the solution X is stored in C. 

COMMON blocks ; None 

Error message: If the reduction to Schur form fails, the message "IN BARSTW, 

EITHER THE SUBROUTINE AXPXB OR ATXPXA WAS UNABLE TO REDUCE A OR B TO 
SCHUR FORM" is printed, and the program is returned to the calling point. 

Field length : 521 octal words (337 decimal) 

Subroutines employed by BARSTW: PRNT, TRANP, EQUATE, NULL, JUXTR, AXPXB, 

ATXPXA, LNCNT 

Subroutines employing BARSTW : VARANCE, CSTAB, DSTAB, RICTNWT, EXPMDFL 

Comments : If a set of equations is to be solved for a collection of 

C matrices with the same A and B or the Schur forms are desired along 
with X, the subroutines AXPXB and ATXPXA should be applied directly. 
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Subroutine TESTSTA 


Description ; The purpose of TESTSTA is to compute and test the eigenvalues 
of a constant real matrix A for stability relative to a parameter a in 
either the continuous or digital sense. In the continuous case, the 
matrix A is classified as stable if the real part of each eigenvalue is 
strictly less than a. Otherwise, A is classified as unstable relative 
to a. In the discrete case, the matrix A is classified as stable if the 
complex modulus of each eigenvalue is less than ct. Otherwise, A is 
declared unstable relative to a. 

Source of software ; Ernest S. Armstrong, LaRC 

Calling sequence; CALL TESTSTA ( A, NA, ALPHA, DISC, STABLE,IOP, DUMMY) 


Input arguments ; 

A Square real matrix packed by columns in a one-dimensional array; 

not destroyed upon return 

NA Two-dimensional vector giving the number of rows and columns of A; 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

ALPHA Scalar a used for testing relative stability 

DISC Logical variable; 

TRUE Test for stability in the discrete sense. 

FALSE Test for stability in the continuous sense. 

lOP Print parameter ; 

0 Do not print results. 

Otherwise Print A, the eigenvalues of A, and stability 
results. 

DUMMY Vector of working space for computations with dimension at least 
(n^ + 5n) where n is the order of A 

Output arguments ; 

STABLE Logical variable; upon normal return, STABLE = TRUE if the 
stability tests are satisfied; otherwise, STABLE = FALSE. 

DUMMY Upon a normal return, the eigenvalues of A are stored in the 
first 2n elements of DUMMY and packed by columns as a 
n X 2 matrix. For DISC = TRUE, the moduli of the eigenvalues 
of A appear in the next n elements. 

COMMON blocks: None 
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Error message; If the computation of the eigenvalues of A fails, the mes- 
sage "IN TESTSTA, THE EIGENVALUE OF A HAS NOT BEEN FOUND AFTER 

30 ITERATIONS" is printed, and the program is returned to the calling 
point. 

Field length ! 366 octal words (246 decimal) 

Subroutines employed by TESTSTA ; EQUATE, EIGEN, JUXTC, PRNT 
Subroutine employing TESTSTA ; ASYMREG 
Comments; None 
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Subroutine EXPSER 


Description ; The purpose of EXPSER is to evaluate the matrix exponential 
eAT 

for a real square matrix A and scalar T. Computation is based on the 
finite-series algorithm described by Kallstrom (ref. 20). The matrix AT 
is scaled by 1/2^ where k is a positive integer chosen so that the scaled 
matrix has eigenvalues within the unit circle in the complex plane. The 
series algorithm is applied to the scaled matrix until the series converges. 
Numerically, convergence is assumed to have occurred when the improvement 
in the element of the finite series of largest magnitude (measured relatively 
if the magnitude is less than unity, and absolutely otherwise) is less than 
the parameter SERCV found in the COI^ON block CONV of subroutine RDTITL. 
Finally, the desired matrix exponential is reconstructed from the exponential 
of the scaled matrix. 

Source of software ; Ernest S. Armstrong, LaRC 

Calling sequence; CALL EXPSER ( A, NA, EXP A, NEXPA,T,I0P, DUMMY) 


Input arguments ; 

A Square matrix packed by columns in a one-dimensional array; not 

destroyed upon return 

NA Two-dimensional vector giving the number of rows and columns of A: 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

T Scalar parameter 

lOP Print parameter; 

0 Do not print results. 

Otherwise Print A, T, and e*"^. 

DUMMY Vector of working space for computations with dimension at 

least 2n^ where n is the order of A 

Output arguments ; 

EXPA Upon a normal return, a square matrix packed by columns in a 

one-dimensional array containing the matrix exponential e*"^ 

NEXPA Two-dimensional column vector giving the number of rows and columns 
of e^"^: 

NEXPA(I) = NA(1) 

NEXPA(2) = NA(2) 
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COMMON block: CONV 


Error megsages : 

(1) The integer k is tested, and if found to be negative, the message 

"ERROR IN EXPSER, K IS NEGATIVE" is printed, and the program returned 
to the calling point. 

(2) If k increases to 1000, the message "ERROR IN EXPINT, K = 1000" is 

printed, and the program is returned to the calling point. 

Field length ; 673 octal words (443 decimal) 

Subroutines employed by EXPSER: MAXEL, UNITY, TRCE, EQUATE, NORMS, SCALE, 

ADD, MULT, PRNT 

Subroutines employing EXPSER: TRANSIT, SAMPL, CNTNREG 

Comments: None 
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Subroutine EXPADE 


Description ; The purpose of EXPADE is to compute the matrix exponential e^, 
where A is a real square matrix stored as a variable-dimensioned two- 
dimensional array. Computation is by the method of Pade approximation 
(ref. 21). The matrix is first scaled by a power of 2 chosen so that the 
eigenvalues of the scaled matrix are within the unit circle in the complex 
plane. The exponential is computed for this scaled matrix using the approxi- 
mation given by the ninth diagonal term in the Pade table for exponential 
approximations. The exponential for the original matrix is then recon- 
structed from the exponential of the scaled matrix. 

Sou rce of software ; LaRC Analysis and Computation Division subprogram library 
with modifications by Ernest S. Armstrong, LaRC 

Calling sequence; CALL EXPADE(MAX,N,A,EA,IDIG,WK,IERR) 


Input arguments; 


MAX 

Maximum first dimension of A as given in the DIMENSION statement 
of the calling program 

N 

Order of matrix A 


A 

Matrix stored in a two-dimensional array with first 
and second at least N; not destroyed upon return 

dimension MAX 

IDIG 

An estimate of the number of accurate digits in the 
ments of absolute value of e^ 

largest ele- 

WK 

Vector of working space for computations, dimensioned at least 
n^ + 8n where n is the order of A 

Output 

arguments: 


EA 

Matrix stored in a two-dimensional array with first 
and second at least N. Upon normal return, EA 
matrix exponential e^ . 

dimension MAX 
contains the 

lERR 

Error parameter; 

0 Normal return 



1 The sum of the absolute values of the elements of A is 

too large for any accuracy, 

2 The Pade denominator matrix is singular. Singularly should 

not occur with exact arithmetic. 

COMMON blocks ; None 

Error messages ; None. The user should examine the parameter lERR upon 
return. 
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Field length; 525 octal words (3**1 decimal) 


Subroutine employed by EXPADE ; 
Subroutines employing EXP APE ; 
Comments: None 


GAUSEL 

None 
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Subroutine EXPINT 

Description ; The purpose of EXPINT is to compute both the matrix exponential 
eAT 

and the integral 

eAs ds 

for a square real matrix A and scalar T. Computation is based on the 
finite-series algorithm described by Kallstrom (ref. 20). The matrix AT 
is scaled by 1/2^ where k is a positive integer chosen so that the 
scaled matrix has eigenvalues within the unit circle in the complex plane. 

The series algorithms are applied to the scaled matrix until convergence 
occurs. Numerically, convergence is assumed to have occurred in each series 
when the improvement in the element of the finite series of largest magnitude 
(measured relatively if the magnitude is less than unity, and absolutely 
otherwise) is less than the parameter SERCV found in the COMMON block CONV 
of subroutine RDTITL. Finally, the desired matrix exponential and the 
integral are reconstructed from the results using the scaled matrix. 

Source of software: Ernest S. Armstrong, LaRC; based on a similar routine in 

“VA'SP (ref. 16) 

Calling sequence; CALL EXPINT(A,NA,B,NB,C,NC,T,IOP, DUMMY) 


Input arguments; 


A Square matrix packed by columns in a one-dimensional array; not 

destroyed upon return 

NA Two-dimensional vector giving the number of rows and columns of A: 

NA( 1 ) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

T Scalar parameter 


lOP 


DUMMY 


Print parameter: 

0 Do not print results. 

Otherwise Print A, T, eAT, g^d 



eA® ds 


Vector of working space for computations, dimensioned at least 2n^ 
where n is the order of A 
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Output arguments; 


B Upon normal return, square matrix packed by columns in a one- 

dimensional array containing the matrix exponential 

eAT 

NB Two-dimensional column vector containing the row and column size 

of e*^: 

NB(1) = NA(1) 

NB(2) = NA(2) 

C Upon normal return, square matrix packed by columns in a one- 

dimensional array containing the matrix 

) e*® ds 

NC Two-dimensional column vector giving the number of rows and columns 

of 


r e*® ds 
*^0 


That is, 

NC(1) = NA(1) 

NC(2) = NA(2) 

COMMON block ; CONV 

Error messages ; 

(1) If k is found to be negative, the message "ERROR IN EXPINT, K IS 

NEGATIVE" is printed, and the program is returned to the calling point. 

(2) If k increases to 1000, the message "ERROR IN EXPINT, K = 1000" is 
^ printed, and the program is returned to the calling point. 

Field length ; 744 octal words (484 decimal) 

Subroutines employed by EXPINT: NORMS, SCALE, UNITY, ADD, EQUATE, MULT, LNCNT, 

PRNT, MAXEL 

Subroutines employing EXPINT; TRANSIT, SAMPL 


Comments: None 
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Subroutine VARANCE 


Description ; The purpose of VARANCE is to compute the steady-state variance 
matrix of the state of the continuous or discrete linear time- invariant 
system; 

Continuous 


x(t) = A x(t) + G Ti(t) 

Discrete 

x(i + 1) = A x(i) + G n(i) 


where A is an n x n asymptotically stable matrix, G is an n x m (m ^ n) 
matrix, and n is a zero-mean white-noise process with continuous intensity 
or discrete variance Q. Following reference 4, the intensity of a continuous 
white-noise process is defined as the coefficient matrix of the 6-function 
in the covariance formula. The steady-state variance matrices, each denoted 
by W, satisfy the Liapunov equations; 

Continuous 


AW + WA* = -GQG' 


Discrete 


W = AWA* + GQG* 


The program provides the option of solving the continuous Liapunov equation 
by either subroutine BILIN or subroutine BARSTW. The discrete Likunov equa- 
tion is solved by using subroutine SUM. The computational parameter EPSA 
used by BARSTW is specified internally as EPSAM found in the COMMON block TOL 
of subroutine RDTITL, 

Source of software ; Ernest S. Armstrong, LaRC 

Calling sequence; CALL VARANCE(A,NA,G,NG,Q,NQ,W,NW,IDENT, DISC, lOP, DUMMY) 


Input arguments ; 


A, G, Q Matrices packed by columns in one-dimensional arrays; not 

destroyed upon normal return except for Q which is replaced 
by GQG*. Storage for Q should be prescribed accordingly in 
the calling program. 
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NA, NG, NQ 


Two-dimensional vectors giving number of rows and columns of 
the respective matrices; for example, 

NA(1) □ Number of rows' of A 
NA(2) = Number of columns of A 
Not destroyed, upon return 

IDENT Logical variable: 

TRUE If G is an identity matrix 

FALSE Otherwise ■ ' 

If IDENT = TRUE, no data are required in G and NG, but 
the variables must still appear as arguments of the calling 
sequence . 

DISC Logical variable: 

TRUE If the discrete case is solved 
FALSE For the continuous case 

lOP Three-dimensional option vector: 

I0P(1) =0 Do not print results. 

Otherwise Print A, G, GQG', Q, and W. 

I0P(2) =0 Do not print from Liapunov equation sub- 

routines employed. 

Otherwise Print from these subroutines. 

I0P(3) = 0 and Solve the Liapunov equation by subroutine 
DISC = FALSE BARSTW. 

I0P(3) ^ 0 and Solve the Liapunov equation by subroutine 
DISC = FALSE BILIN. 

I0P(3) is not required if DISC = TRUE. 

DUMMlf Vector of working spaces for computations dimensioned at least: 


2(n 

+ 1)2 

for 

DISC = FALSE 

and 

I0P(3) = 0 

i»n2 

+ 2n 

for 

DISC = FALSE 

and 

I0P(3) ^ 0 

4n2 


for 

DISC = TRUE 




Output arguments : 

W Upon normal return, matrix packed by columns in a one- 

dimensional array holding the steady-state variance matrix 

NW Two-dimehsional vector giving number of rows and columns of W: 

upon normal return, 

NW(1) = NA(1) 

NW(2) = NA(2) 

COMMON block: TOL 
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Error messages; None 


Field length ; 513 octal words (331 decimal) 

Subroutines employed by VARANCE ; PENT, LNCNT, MULT, TRANP, BILIN, BARSTW, SUM 
Subroutines employing VARANCE ; None 

Comments ; In determining the storage allocation requirements for VARANCE, it 
was assumed that m S n. For m > n, the matrix GQG' should be computed 
externally and input as Q with IDENT = TRUE. 
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Subroutine CTROL 


Description; The purpose of CTROL is to evaluate the controllability matrix 


C = [b, AB, . . A^-Ib] 


for a real constant (A,B) pair. The matrix A is n x n and B is 
n X r with r = n. Options are provided to compute both the rank and singu- 
lar values of C along with the controllability canonical form (ref. 4) for 
the (A,B) pair. For the optional computations, the singular-value decom- 
position algorithm found in subroutine SNVDEC is applied to factor 


C' = [b, AB, . . ., An-i'B] 
as 

c' = UQV 
or 

C = VQU' 


The number of nonzero elements on the diagonal of the n x n matrix Q 
determines the rank of C. Assuming C has rank = n, the first ^ col- 
umns of V form an orthonormal basis for the range space of C and the next 
(S' + 1 ) to n columns form an orthonormal basis for the orthogonal complement 
to the range space of C. Hence, the pair (V*AV,V*B) represents the con- 
trollability canonical form for the original (A,B) pair. 

Source of software ; Ernest S. Armstrong, LaRC 

Calling sequence; CALL CTROL(A,NA,B,NB,C,NC, TOP, IAC,IRANK, DUMMY) 


Input arguments ; 

A, B Matrices packed by columns in one-dimensional arrays; not destroyed 

upon return 

NA, NB Two-dimensional vectors giving the number of rows and columns of 
respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 
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lOP Five-dimensional option vector: 


lOP(l) = 0 
otherwise 

Do not print A, B, and C. 
Print A, B, and C. 

I0P(2) = 0 
Otherwise 

Return after computing C. 
Compute rank of C. 

I0P(3) = 0 
Otherwise 

Do not print rank of C, zero test employed to 
determine rank, and singular values of C. 
Print these data. 

I0P(4) = 0 
Otherwise 

Return after rank computation. 

Compute the controllability canonical form. 

I0P(5) = 0 
Otherwise 

Return without printing controllability form. 
Print V'AV, V'B, and V before returning. 


lAC Parameter used to specify zero test for rank computation. 

Singular values are considered zero if they do not exceed 

ZTEST = lie II X using the matrix Euclidean norm. 

DUMMY Vector of working space for computations dimensioned at least: 
nr(n+1) = K To compute C only 

max [K,n2+n+nr(n-r+1 )3 □ L To compute rank of C 
max(L,3n^) To compute controllability 

canonical form 


Output arguments : 

C Matrix packed by columns in a one-dimensional array dimensioned 

at least n^r. Upon normal return, C contains the controlla- 
bility matrix for the (A,B) pair. 

NC Two-dimensional vector giving the number of rows and columns of C: 

upon normal return, 

NC(1) = NA(1) = n 

NC(2) = NA(1) X NB(1) = nr 

IRANK Upon normal return, scalar giving the rank of C 

DUMMY Upon normal return for rank computation only, the first nr(n-r+1) 

elements contain the matrix U, the next n elements contain the 
singular values of c’ , and the next n^ elements contain the 
matrix V. Upon normal return for the canonical form computation, 
the first n^ elements contain the matrix V’, the next nr ele- 
ments contain the matrix V’B, and, after the first 2n^ elements 
of DUMMY, the next n^ elements contain the matrix V*AV. All 
matrices are packed by columns into the corresponding sections 
of the one-dimensional array DUMMY. 
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COMMON blocks ; None 
Error messages ; 

(1) If SNVDEC fails to compute the singular values of C, the message "IN 

CTROL, SNDVEC HAS FAILED TO CONVERGE TO THE SINGULAR VALUE AFTER 

30 ITERATIONS" is printed, and the program is returned to the calling 
point. 

(2) If an eigenvalue of Q is greater than but close to ZTEST (see sub- 

routine SNVDEC), the message "IN CTROL, THE MATRIX SUBMITTED TO SNVDEC 
USING ZTEST = IS CLOSE TO A MATRIX WHICH IS OF LOWER RANK / IF 

THE ACCURACY IS REDUCED THE RANK MAY ALSO BE REDUCED/CURRENT 
RANK = " is printed, and the computation continues. 

Field length ; 656 octal words (430 decimal) 

Subroutines employed by CTROL ; EQUATE, MULT, JUXTC, PRNT, LNCNT, TRANP, SNVDEC 
Subroutines employing CTROL ; None 

Comments ; Employing duality theory, CTROL may also be used to compute the 
reconstructibility canonical form (ref. 4). By combining controllability 
and reconstructibility results from CTROL, the full canonical structure 
theorem (ref. 22) can be implemented. 
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Subroutine TRANSIT 


Description ; The purpose of TRANSIT is to compute and print the transient 
response of the time-invariant continuous or discrete system. 

Continuous 


X(t) = AX(t) + B U(t) 


Discrete 


X(i + 1) = AX(i) + B U(i) 
together with Y and U where 
Y = HX + GU 
U = -FX + V 


for given 


X(0) = Xo 


from the initial time or stage zero to an input final time or stage. 

The matrices A, B, X, and H are dimensioned n x n, n x r (r $ n), 
n X p (p g n), and m x n (m S n) , respectively, with the other matrices com- 
patible. If the matrix (A - BF) is asymptotically stable in the appropriate 
continuous or discrete sense, the steady-state value of X given by 

Continuous 


X = -(A - BF)-"'bv 


Discrete 


X = [l - (A - BF)]“^BV 

where I is an identity matrix, is computed and printed except when V is 
a null matrix. 
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Source of software; Ernest S. Armstrong, LaRC 


Calling sequence; CALL TRANSIT(A,NA,B,NB,H,NH,G,NG,F,NF,V,NV,T,X,NX,DISC, 
STABLE, lOP, DUMMY) 


Input arguments ; 

A, B, H, Compatible matrices packed by columns in one-dimensional 

G, F, V arrays; not destroyed upon normal return 


NA, NB, NH, Two-dimensional vectors giving the number of rows and columns 
NG, NF, NV of respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 


T Two-dimensional vector; 

T(1) is the final time in the continuous case given as an 
integer multiple of T(2). 

T(2) is a time point within [0,T(1)] ; printing is at 
multiples of T(2). 

The input T is not required if the response to a discrete 
system is to be computed, but the argument must still 
appear in the calling sequence; not destroyed upon normal 
return. 


X 


NX 


DISC 


STABLE 


Matrix packed by columns in a one-dimensional array contain- 
ing the initial value Xg at time or stage zero. X is 
destroyed upon normal return. 

Two-dimensional vector giving the number of rows and columns 
of X; 

NX(1) = NA(1) 

NX(2) = NU(2) 

Not destroyed upon normal return 
Logical variable; 

TRUE If the response to the discrete system is required 
FALSE For the response to the continuous system 

Logical variable; 

TRUE If (A - BF) is asymptotically stable in 

the appropriate sense 
STABLE = FALSE Otherwise 
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lOP 


Four-dimensional 

option vector; 

lOP(l) 

= 0 


If 

H is i 

lOP(l) 

= 1 


If 

H is j 

lOP(l) 

# 0 

or 1 

For 

other 

I0P(2) 

= 0 


If 

G is ; 

IOP(2) 

= 1 


If 

G is i 

I0P(2) 

4 0 

or 1 

For 

other 

I0P(3) 

= 0 


If 

V is i 

I0P(3) 

= 1 


If 

V is ; 

IOP(3) 

* 0 

or 1 

For 

other 

IOP(4) 

is 

the 

terminal stage 


H matrices 


a null matrix 


G matrices 


V matrices 


r the response to a discrete 
system; not ^required if the continuous system response is 
computed (DISC = FALSE) . 


DUMMY 


Vector of working space for computations, dimensioned at 
least Tn^ where n is the order of A. 


Output arguments; 


X 


Upon normal return, the value of X at time T(1) or 
state I0P(4) 


DUMMY Upon normal return, the first np elements of DUMMY contain 

the steady-state value of X, packed by columns in a one- 
dimensional array, when STABLE = TRUE and I0P(3) 4 0« 

COMMON blocks ; None 

Error message ; When computation of the continuous steady-state X values is 
attempted and (A - BF) is found to be singular, the message "IN TRANSIT, 
THE MATRIX A-BF SUBMITTED TO GELIM IS SINGULAR" is printed, and the program 
is returned to the calling point. For the discrete case, a similar message 
concerning the matrix p - (A - BF)] is also printed. 

Field length ; 1530 octal words (856 decimal) 

Subroutines employed by TRANSIT; PRNT, EXPSER, EXPINT, MULT, EQUATE, LNCNT, 
SCALE, ADD, TRANP, GELIM, UNITY, SUBT 


Subroutines employing TRANSIT ; None 

Comments ; When the matrices H, G, or V take on their special null or 
identity matrix values, no data need be input, but they must still appear as 
arguments of the calling sequence. 
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PBIMARY SUBROUTINES FOR IMPLEMENTING LQG CONTROL LAW DESIGN 


Subroutine SAMPL 

Description; The purpose of SAMPL is to compute the matrix functions, 


Q(T) 


= r eA''^Qe^'^ dT 

W(T) = 2 f eA'o H(T,0) dx 
^0 

= f [r + H'(T,0) Q H(T,0)] dX 


R(T) 


where 


.0, = r 


H(t,0) = \ eA'^ B dX 
^0 


for constant real matrices A, B, Q=Q'S0, R'=R>0, and scalar 

T > 0. The dimensions of A and B are n x n and n x r (r S n), 
respectively. Other matrices have compatible dimensions. The^program has 
the option of computing Q(T) without evaluating W(T) and R(T). 

The matrix Q is the reconstructibility Gramian (ref. 4) for the system 
x(t) = A x(t) 

y(t) = D x(t) > (0 ^ t ^ T) 

Q = D'D ^ 

The set of matrices Q, W, and R occur naturally in the optimal sampled- 
data linear regulator problem (refs. 23 and 24). Given the time-invariant 
linear system, 


x(t) = A x(t) + B u(t) 
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the optimal sampled-data regulator problem (OSR) occurs when u(t) 
(0 S t S tf < °°) is required to minimize 


(tf) S x(tf) + r [x'(t) q x(t ) + u'(t) r u(t)]| dr 


subject to the restrictions that u(t) be constant over subintervals 
[ti.ti+l] (i = 0,1,...,N-1), 0< to < ti < t2 . . . < tjj = tf within 
the interval []0,tfj, and S = S' ^ 0. The OSR problem transforms 
directly into an optimal discrete regulator problem in which uCt^) 

(i = 0,1,..., N-1) is chosen to minimize 


J = x'(tf) S x(tf) 


N-1 

It 


x'(ti) Q(Ati) x(ti) 


i=0 


+ x'(ti) VKAtf) u(ti) + u'(ti) R(Ati) u(ti)J 


with 


Ati = tj^+1 - tj^ 


and 


x(ti+i) = e^^i x(ti) + H(Ati,0) uCt^) 


The subroutine SAMPL computes the weighting matrices for the OSR problem for 
the case At^ = T. 

Computation is based on the finite-series algorithm described by Armstrong 
and Caglayan (ref. 13). The matrix AT is scaled by 1/2^ where k is a 
positive integer chosen so that the scaled matrix has eigenvalues within the 
unit circle in the complex plane. The algorithm (ref. 13) is applied to the 
scaled matrix until convergence occurs. Numerically, convergence is assumed 
to have occurred in each series when the improvement in the element of largest 
magnitude (measured relatively if the magnitude is less than unity, and abso- 
lutely otherwise) is less than the parameter SERCV found in the COMMON block 
CONV of subroutine RDTITL. Afterward, the desired OSR weighting matrices are 
reconstructed from the results with the scaled AT matrix, as described in 
reference 13. 

Source of software ; Ernest S. Armstrong, LaRC 

Calling sequence: CALL SAMPL(A, NA,B,NB,Q,NQ,R, NR, W,NW,T,I0P, DUMMY) 
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Input arguments; 


A, B, Q, R Matrices packed by columns in one-dimensional arrays. A 

and B are not^destroyed upon return, but Q and R 
are. If only Q is to be computed, data for B and R 
need not be input, but the related arguments should still 
appear in the calling sequence. 

NA, NB, NQ, NR Two-dimensional vectors giving the number of rows and 

columns of respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon normal return 

T Positive scalar parameter 

lOP Two-dimensional option vector: 

I0P(1) =0 Do not print results of computation. 
Otherwise Print input and computed results. 

I0P(2) = 0 Solve for the reconstructibility Gramian 
only. 

Otherwise Solve for the complete set of OSR weighting 
matrices. 

DUMMY Vector of working space for computations with dimension 

at least: 

4n2 for I0P(2) r 0 

7n2 for I0P(2) f 0 

Output arguments : 

Q Upon normal return, the matrix Q(T) is stored in Q and 

packed by columns in the one-dimensional array. 

W Matrix packed by columns in a one-dimensional array with 

dimension at least nr. Upon normal return, W contains 
the matrix W(T) if computation is required. If the 
computation of W(T) is not required, W is not needed, 
but the argument should still appear in the calling 
sequence . 

NW Two-dimensional vector giving, upon normal return, the 

number of rows and columns of W(T): 

NW(1) = NA(1) 

NW(2) = NB(2) 

R Upon normal return, the matrix R(T) is stored in R and 

packed by columns in the one-dimensional array. 
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COMMON block ; CONV 

Error message ; If k reaches 1000 in the scaling of AT, the message "ERROR 
IN SAMPL, K = 1000" is printed, and the program is returned to the calling 
point . 

Field length ; 2521 octal words (1361 decimal) 

Subroutines employed by SAMPL: PRNT, LNCNT, NORMS, SCALE, EQUATE, MULT, TRANP, 

ADD, MAXEL, EXPSER, EXPINT 

Subroutines employing SAMPL : None 

Comments : The variance matrix for the state of the time-invariant linear 

system 

x(t) = A x(t) + B U(t) 


driven by white noise U with intensity V is given at time T by 


G(T) = eAT G(0) e^'T + J* eA^^BVE' e^ dT 


and can be computed through the matrix exponential and SAMPL subroutines of 
ORACLS. For the second term in G(T), use SAMPL and compute the reconstructi- 
bility Gramian with A replaced by A’ and Q □ BVB' . 
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Subroutine PREFIL 


Description : The purpose of PREFIL is to compute an r x n (r ^ n) matrix F 

which, when used in the vector equation, 


u = -Fx + V 


eliminates the cross-product term in the quadratic scalar function, 


x^Qx + x’Wu + u’Ru 


where Q = Q’ ^ 0, W, and R = R* > 0 are constant matrices. Specifically, 
1 ^ 

F = R“^ - 
2 


and, after substitution, the quadratic function becomes 


x’Qx + v*Rv 


where 


Q n Q - 


W 

- F 
2 


If the transformation is also applied to a continuous or discrete linear 
time-invariant system, 

Continuous 


X = Ax + Bu 


Discrete 


x(i+1) = A x(i) + B u(i) 


the closed-loop response matrix is 


A = A - BF 


Options are provided within PREFIL to compute both Q and A in addition 
to F. 
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Source of software ; Ernest S. Armstrong, LaRC 

Calling sequence ; CALL PREFIL(A,NA,B,NB,Q,NQ,W,NW,R , NR ,F ,NF ,IOP , DUMMY) 

Input arguments ; 

A, B, Q, Matrices packed by columns in one-dimensional arrays. 

W, R Inputs B and R are not destroyed upon normal return. 

If only the matrix F is required, input A, B, and Q 
are not required but must still appear in the calling 
sequence. Similarly, if only F and Q are required, 
input B is not required, but B must still appear in the 
calling sequence. 

NA, NB, NQ, Two-dimensional vectors giving the number of rows and columns 

NW, NR in the respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon return 

lOP Three-dimensional option vector: 

I0P(1) =0 Do not print results. . 

Otherwise Print input data, F, and (q, a) when 

computed. 

I0P(2) =0 Do not compute 

Otherwise Compute F and Q and return. 

I0P(3) H 0 Do not compute A^ 

Otherwise Compute F and A and return. 

DUMMY Vector of working space for computations, dimensioned at 

least n^ + r where n and r are the order of A and R, 
respectively 

Output arguments : 

A If A is computed, the A array contains, upon normal return, 

the matrix A packed by columns. 

Q If Q is computed, the Q array contains, upon normal return, 

the matrix Q packed by columns. 

F Upon normal return, the matrix F packed by columns in a one- 

dimensional array of dimension at least nr 

NF Upon normal return, a two-dimensional vector giving the number 

of rows and columns of F: 

NF(1) = NR(1) 

NF(2) = NA(1) 
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COMMON blocks: None 


Error message: If the matrix R is found not to be symmetric positive definite 

by the subroutine SYMPDS, the message "IN PREFIL, THE MATRIX R IS NOT 
SYMMETRIC POSITIVE DEFINITE" is printed, and the program is returned to the 
calling point. 

Field length : 457 octal words (303 decimal) 

Subroutines employed by PREFIL: PRNT, LNCNT, TRANP, SCALE, EQUATE, SYMPDS, 

MULT, SUBT 

Subroutine employing PREFIL : IMPMDFL 

Comments : Subroutines which follow provide solution algorithms for LQG problems 

having a quadratic performance index without cross-product terms. Problems 
with such terms may be transformed into equivalent problems without cross 
products by performing a control variable transformation (ref. 4): 


u = -Fx + V 


where u is the original control, x is the state vector, and v is the 
new control with F computed from PREFIL. 



Subroutine CSTAB 


Description ; The purpose of CSTAB is to compute a gain matrix F which, when 
used in the control law 

u = -Fx 

and applied to the stabilizable linear time- invariant continuous system, 

X = Ax + Bu ' 


produces a closed-loop response matrix (A - BF) whose eigenvalues lie in 
the complex left half-plane. The primary use for CSTAB in ORACLS is to gen- 
erate a stabilizing gain matrix for initializing the quasilinearization 
method for solving the continuous steady-state Riccati equation (ref. 9). 

The matrix (A - BF) computed here has some interesting root- locus proper- 
ties (ref. 25) which may make the control law applicable in other areas. 

Computation follows the method described by Armstrong (ref. 11). The 
matrices A and B are of dimension n x n and n x r (r ^ n), respec- 
tively. A scalar parameter 3 > 0 is first selected so that the matrix 


A = -(A + 31) 


has eigenvalues in the complex left half-plane, where I is the identity 
matrix. The matrix is used to form a Liapunov equation 


AZ + ZA' = -2BB* 

whose solution Z is used to compute F through 


F = B'Z+ 


where + denotes matrix pseudoinverse. It can be shown (ref. 25) that 
Re(xA-BF) _ .3 


where X^'BF controllable eigenvalue of the (A,B) system. The 

option is provided to input 3 directly or to have 3 computed internally 
as 
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B = s 



\ 

1 f Pi 
max Re\Xj^ 

_ i 

j + 0.001 


A 

where s is an input scale factor and 
the eigenvalue computation in CSTAB fails, 


(i " 1|2,«*«|n) 


are the eigenvalues of A. 
6 is set to 


If 


6 = 2||All 


using the it-) matrix norm (see subroutine NORMS). The Liapunov equation 
for Z can be solved using either subroutine BILIN or subroutine BARSTW. 
The pseudoinverse of Z is computed through subroutine SNVDEC. Computa- 
tional parameters for BILIN, BARSTW, and SNVDEC are provided internally 
through the COMMON block TOL found in subroutine RDTITL. 

Source of software ; Ernest S. Armstrong, LaRC 

Calling sequence; CALL CSTAB(A,NA,B,NB,F,NF,IOP,SCLE, DUMMY) 


Input arguments ; 

A, B Matrices packed by columns in one-dimensional arrays; not destroyed 

' upon normal return 


NA, NB Two-dimensional vectors giving the number of rows and columns of 
the respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon normal return 


lOP 


SCLE 

DUMMY 


Three-dimensional option vector: 

I0P(1) =0 Do not print results. 

Otherwise Print input, 3, F, and eigenvalues of (A-BF). 


I0P(2) =0 Do not compute parameter 3 but use 3 = SCLE. 
Otherwise Compute 3 using s = SCLE. 


I0P(3) = 0 Use the BARSTW algorithm to solve the Z equation. 
Otherwise Use BILIN. 


Parameter used to define 3 

Vector of working space for computations of dimension at least: 
2n2 + 2(n +1)* if I0P(3) = 0 
6n2 + 2n if I0P(3) ^ 0 
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Output arguments; 


F Matrix packed by columns into a one-dimensional array of dimension 

at least nr. Upon normal return, F contains B'Z+. 

NF Two-dimensional vector holding, upon normal return, the number of 

rows and columns of F: 

NF(1) □ NB(2) 

NF(2) = NA(1) 

COMMON block ; TOL 

Error messages ; 

(1) In the 3 computation section, if the eigenvalue computation for A 

fails, the message "IN CSTAB, THE SUBROUTINE EIGEN FAILED TO DETERMINE 

THE EIGENVALUE FOR THE MATRIX A AFTER 30 ITERATIONS" is 

printed, and the computation continues with 3 = 2||A||. 

(2) If SNVDEC fails to compute the singular values of Z, the message "IN 

CSTAB, SNVDEC HAS FAILED TO CONVERGE TO THE SINGULAR VALUE 

AFTER 30 ITERATIONS" is printed, and the program is returned to the 
calling point. 

(3) If a singular value of Z is greater than but close to the ZTEST value 

used (see SNVDEC and TOL), the message "IN CSTAB, THE MATRIX SUBMITTED 
TO SNVDEC USING ZTEST = IS CLOSE TO A MATRIX OF LOWER RANK/IF 

THE ACCURACY lAC IS REDUCED THE RANK MAY ALSO BE REDUCED/CURRENT 

RANK = " is printed along with the singular values of Z after 

which computation continues. 

Field length ; 1076 octal words (574 decimal) 

Subroutines employed by CSTAB: EQUATE, EIGEN, LNCNT, TRANP, MULT, SCALE, 

BARSTW, BILIN, SNVDEC, PRNT, SUBT, NORMS 

Subroutines employing CSTAB; DSTAB, ASYMREG 


Comments: None 
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Subroutine DSTAB 


Description ; The purpose of DSTAB is to compute a gain matrix F which, when 
used in the control law, 

u = -Fx 


and applied to the stabilizable linear time- invariant discrete system. 


x(i+1) = A x(i) + B u(i) 


produces a closed-loop response matrix (A - BF) whose eigenvalues lie inside 
the unit circle of the complex plane. The primary use for DSTAB in ORACLS 
is to generate a stabilizing gain matrix for initializing the quasilineariza- 
tion method for solving the discrete steady-state Riccati equation (ref. 10). 
The matrix (A - BF) computed here has some interesting root-locus properties 
(ref. 25) which may make the control law applicable in other areas. 

Computation follows the method described by Armstrong and Rublein (ref. 12). 
The matrices A and B are of dimension n x n and n x r (r = n), 
respectively. A scalar a is first selected so that 


0 < a < min ( 1 , min 



(i H 1,2, ...,n) 


where | | denotes the nonzero complex moduli of eigenvalues of A. Next, 
solve the matrix equation, 


AZA’ = a2z + 2BB’ 


for Z = Z' i 0, and by assuming the controllable eigenvalues of A are 
nonzero, the stabilizing gain is given by 


F = B'(Z + BB’)+A 


where + denotes matrix pseudo inverse. If any controllable eigenvalue 
of A is zero or unknown, it can be made nonzero by applying to the system 
a prefilter with gain G so that the controllable eigenvalues of (A - BG) 
are strictly in the complex left half-plane. Such a gain can be computed 
from CSTAB. Afterward, DSTAB is applied to the (A - BG,B) pair to compute 
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a digital stabilizing gain F with the net stabilizing gain F for the 
original system given by 


F = G + F 


The matrix equation in Z is transformed into 


AZ + ZA' = BB' 


with 


A = (al + A)-''(A - al) 
B = 2(al + A)-''b 


and solved by subroutine BARSTW. The symbol I denotes an identity matrix. 

The option is provided to input a directly or to have a computed 
internally as 


a = s min 


min 

i 



(i = 1,2 n) 


A-BG 

where are nonzero complex moduli of eigenvalues of (A - BG), 

0 < s < 1 is input scale factor, and G is a gain computed, if requested, 
internally from CSTAB to cause the controllable eigenvalues of (A-BG) to 
be nonsingular. The pseudoinverse is computed through SNVDEC. Computational 
parameters for SNVDEC and BARSTW are provided internally through the COMMON 
block TOL found in subroutine RDTITL. 


Source of software : Ernest S. Armstrong, LaRC 

Calling sequence; CALL DSTAB(A,NA,B,NB,F,NF, SING, IOP,SCLE, DUMMY) 


Input arguments ; 

A, B Matrices packed by columns in one-dimensional arrays; not destroyed 

upon normal return 

NA, NB Two-dimensional vectors giving the number of rows and columns in 
the respective matrices; for example, 

NA(1) B Number of rows of A 
NA(2) = Number of columns of A 
Not destroyed upon normal return 
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SING 


Logical variable : 

FALSE If the controllable eigenvalues of A are known to 
be nonzero : 

TRUE Otherwise, and a gain G will automatically be computed 
(using CSTAB) before proceeding to the digital 
algorithm • 

lOP Two-dimensional option- vector: 

lOP(l) =0 Do not ‘print results. 

Otherwise Print input, F, a, eigenvalues of A - BF, and 
their magnitudes. 

I0P(2) =0 Do not compute parameter a but use a = SOLE. 

Otherwise Compute a by using s = SCLE. 

SCLE Parameter used to define a 

DUMMY Vector of working space for computations, dimensioned at least 

4n^ + 2(n + 1)^ 

Output arguments : 

F Matrix packed by columns into a one-dimensional array of dimension 

at least nr. Upon normal return, F contains B'(Z + BB')'*'A. 

NF Two-dimensional vector holding, upon normal return, the number of 

rows and columns in F: 

NF(1) = NB(2) 

NF(2) = NA(1) 

COMMON block ; TOL 

Error messages : 

(1) In the a computation, if the eigenvalue computation for (A - BG) 

fails, the message "IN DSTAB, THE PROGRAM EIGEN FAILED TO DETERMINE 

THE EIGENVALUE FOR THE MATRIX A - BG AFTER 30 ITERATIONS" is 

printed followed by the matrices (A - BG) and G, if SING = TRUE, 
and the program is returned to the calling point. 

(2) If the matrix A + al (or (A - BG) + al) is found to be singular, 

the message "IN DSTAB, GELIM HAS FOUND THE MATRIX A + (ALPHA)I 
(or (A-BG) + (ALPHA)I) SINGULAR" is printed followed by A, G, 
and a, and the program is returned to the calling point. 

(3) If SNVDEC fails to compute the singular values of (Z + BB*), the 

message "IN DSTAB, SNVDEC HAS FAILED TO CONVERGE TO THE SINGULAR 

VALUE AFTER 30 ITERATIONS" is printed, and the program is returned to 
the calling point. 
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(4) If a singular value of (Z + BB') is greater than, but close to the 

ZTEST value used (see SNVDEC and TOL), the message "IN DSTAB, THE 
MATRIX SUBMITTED TO SNVDEC, USING ZTEST = , IS CLOSE TO A 

MATRIX OF LOWER RANK/IF THE ACCURACY lAC IS REDUCED THE RANK MAY ALSO 

BE REDUCED/CURRENT RANK = " is printed along with the singular 

values of (Z + BB') after which computation continues. 

(5) If the eigenvalue computation for (A - BF) fails, the message "IN DSTAB, 

THE PROGRAM EIGEN FAILED TO DETERMINE THE EIGENVALUE FOR THE 

A-BF MATRIX AFTER 30 ITERATIONS" is printed along with the computed 
eigenvalues, and the program is returned to the calling point. 

Field length ; 1605 octal words (901 decimal) 

Subroutines employed by DSTAB: CSTAB, MULT, SUBT, EQUATE, EIGEN, GELIM, LNCNT, 

PRNT, TRANP, MULT, SCALE, BARSTW, ADD, SNVDEC, JUXTC 

Subroutine employing DSTAB ; ASYMREG 

Comments: None 
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Subroutine DISCREG 


Description ; The purpose of DISCREG is to solve the time-invariant discrete- 
time linear optimal output regulator problem with noise-free measurements. 
Given the digital linear system 


x(i+1) = A x(i) + B u(i) + w(i) 


where x(0) = xq is given, A and B are constant matrices of dimension 
n X n and n x r (r S n), respectively, and w(i) (i = 0,1,..., N-1) is a 
sequence of uncorrelated zero-mean stochastic variables with variance 
matrices V(i). Outputs, or controlled variables, are of the form 


y(i) = H x(i) 


where H is a constant m x n (m ^ n) matrix. The considered optimal regu- 
lator problem is to find the control sequence u(i) which minimizes 


J = E 


N-1 

^ |y'(i+1) Q y(i+1) + u*(i) R u(i)3 
i=0 


x'(N) Pjg x(N) 


with Q = Q' ^ 0, Pjg = Pn' = R = R* > 0. The symbol E denotes 

expected value. The solution is given by (ref. 26) 


u(i) B -F(i) x(i) 1 


F(i) 

P(i) 

<t)(i) 


[r + B' P(i+1) B' P(i+1) A 

ct)'(i) P(i+1) <|)(i) + F'(i) R F(i) + H'QH 
A - BF(i) 


(i = 0,1 


,N-1) 


n 



where P(N) = Pjj. The minimal value of the criterion function is 


N-1 

J = x'(0) P(0) x(0) + 

j=0 



tr [V(j) P(j+1)] 


where tr denotes the trace of a matrix. 

The program DISCREG repetitively evaluates the solution equations from an 
input value N to stage zero or until P(i) converges to a steady-state 
value, whichever occurs first. Numerically, convergence is assumed to have 
occurred when the improvement in the element of largest magnitude (measured 
relatively if the magnitude is less than unity, and absolutely otherwise) is 
less than the parameter RICTCV found in COMMON block CONV of subroutine 
RDTITL. The program DISCREG does not evaluate the optimal value of the 
criterion function or print a response trajectory. The subroutine SNVDEC 
is used to evaluate the matrix inverse in the F(i) equation to allow the 
option of using nonnegative definite (instead of strictly positive definite) 
matrices R. Parameters for SNVDEC are set internally by use of the COMMON 
block TOL of RDTITL. 

Source of software ; Ernest S. Armstrong, LaRC 

Calling sequence: CALL DISCREG(A,NA,B,NB,H,NH,Q,NQ,R,NR,F,NF,P,NP,IOP, 

IDENT, DUMMY) 

Input arguments: 


A, B, H, Matrices packed by columns in one-dimensional arrays. Upon 

Q, R normal return, all except Q are not destroyed. 

NA, NB, NH, Two-dimensional vectors giving the number of rows and columns 
NQ, NR of the respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 

Not destroyed upon normal return except for NQ = NA 

P Symmetric nonnegative definite matrix packed by columns in 

a one-dimensional array. On input, P contains the 
matrix Pjj. Afterward, intermediate values of P(i) are 
stored in the array P. 

NP Two-dimensional vector giving the number of rows and columns 

of P: 

NP(1) = Number of rows in P 
NP(2) = Number of columns in P 
Not destroyed upon normal return 
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lOP 


IDENT 


DUMMX 


Three-dimensional option vector: 

I0P(1) =0 Do not print results. 

Otherwise Print input, and F(i) and P(i) at stage 

zero or steady state, whichever occurs first. 


IOP(2) =0 Do not print at intermediate states between 
initial and final states. 

Otherwise Print stage count and values of F(i) and 
P(i), regardless of printing specified by 
I0P(1). 


I0P(3) is terminal stage N for the optimal regulator 
problem . 

Logical variable: 

TRUE If H is an identity matrix 

FALSE Otherwise 

For IDENT = TRUE, no data need be input for H, but the 
related arguments should still appear in the calling 
sequence . 

Vector of working space for computations of dimension at 
least 4n^ 


Output arguments: 


Q Upon normal return, Q is replaced by H'QH. The array Q 

must be dimensioned at least n^. 

F Matrix packed by columns in a one-dimensional array of dimen- 

sion at least nr. Upon normal return, F contains the 
value of F(i) at stage zero or the stage at which numeri- 
cal steady-state convergence occurs. 

NF Two-dimensional matrix giving, upon normal return, the number 

of rows and columns of F: 

NF(1) = NB(2) 

NF(2) = NA(1) 

P Upon normal return, P contains the value of P(i) at stage 

zero or the stage at which steady-state convergence occurs. 

COMMON blocks : TOL, CONV 

Error messages : 

(1) If SNVDEC fails to compute the singular values of any j^R + B’ P(i) Bj , 

the message ’'IN DISCREG, SNVDEC HAS FAILED TO CONVERGE TO THE 

SINGULAR VALUE AFTER 30 ITERATIONS" is printed, and the program is 
returned to the calling point. 
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(2) If a singular value of jj? + B' P(i) b] is greater than but close to the 
ZTEST value used (see SNVDEC and TOL) , the message "IN DISCREG, THE 
MATRIX SUBMITTED TO SNVDEC USING ZTEST = IS CLOSE TO A MATRIX 

OF LOWER RANK/IF THE ACCURACY lAC IS REDUCED THE RANK MAY ALSO BE 

REDUCED/CURRENT RANK = " is printed along with the singular 

values of fR + B'PCDbJ after which computation continues. 

Field length ; 1355 octal words (7^9 decimal) 

Subroutines employed by DISCREG: LNCNT, PRNT, MULT, TRANP, EQUATE, SNVDEC, 

ADD, SUBT, MAXEL 

Subroutine employing DISCREG : ASYMREG 

Comments : Of course, if w(i) = 0 (i = 0,1,...), the algorithm in DISCREG 

generates the solution to the discrete deterministic optimal linear output 
regulator problem. 
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Subroutine CNTNREG 


Description: The purpose of CNTNREG is to solve the time-invariant continuous- 

time linear optimal output regulator problem with noise-free measurements. 
Given the continuous linear system, 


x(t) = A x(t) + B u(t) + w(t) 


where x(0) = xq is given, A and B are constant matrices of dimension 
n X n and n x r (r S n), respectively, and w(t) is zero-mean Gaussian 
white noise with intensity V(t) . Outputs, or controlled variables, are 
modeled as 


y(t) = H x(t) 


where H is a constant m x m (m S n) matrix. The considered optimal regu- 
lator problem is to find the control function u(t) which minimizes 


J = 



(t) Q y(t) + u'(t) R u(t)] 


dt + x’ (ti ) ?■) x(ti ) 


where Q = Q'S0, P-|=:Pi*^0, and R = R’ > 0. The symbol E denotes 

expected value. The solution is given by reference 4 as 


u(t) = -F(t) x(t) 
F(t) = R“''b' P(t) 


-dP(t) 

dt 


= H’QH + A’ P(t) + P(t) 


A - P(t) BR-1 b' P(t) 


where P(ti) = P-|. The minimal value of the criterion function is 


x'(0) P(0) x(0) + 



tr [P(t) V(t)] dt 


where tr denotes the trace of a matrix. 
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The computational algorithm used in CNTNREG to solve the Riccati equation is 
that due to Vaughan (ref. 8) as described by Kwakernaak and Sivan (ref. 4). 
From 


A 


Z = 


-H'QH 


-BR-^B' 

-A’ 


The Z matrix has the property that if X is an eigenvalue, then -X is 
also an eigenvalue. By assuming that Z has linearly independent eigen- 
vectors, there exists a matrix W such that 


Z = 




w-1 


where A is a diagonal matrix constructed as follows. If an eigenvalue X 
of Z has Re(X) >0, it is a diagonal element of A, and -X is automati- 
cally placed in -A. If Re(X) □ 0, one of the pair (X,-X) is arbitrarily 
assigned to A and the other to -A. The matrix W is composed of eigen- 
vectors of Z; the ith column vector of W is the eigenvector of Z corre- 
sponding to the eigenvalue in the ith diagonal position of diag(A,-A). In 
practice, it is noted that if X is a complex eigenvalue of Z with eigen- 
vector V, then their complex conjugates X and v are also an eigenpair 
and []Re(v) , Im( v)l is placed in W instead of (v,v) in order to avoid 
complex arithmetic. This has the effect of making A block diagonal where 
every (X,X) pair go together to form a 2x2 real matrix placed at the 
former (X,X) entry in A. Next partition the 2n x 2n matrix W into 
four n X n blocks as 



"Wl1 

W 12 

w = 

_W21 

W22_ 




The solution 

P(t) 

of the 


P(t) = [W22 + W 21 G(ti - 

with 

G(t) = 


equation can be written as 
t)] [Wi2 + W 11 G(ti - t^"'' 
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and 


S = -(W21 - PiWii)-''(W 22 - P1W12) 


If A is composed of eigenvalues with strictly positive real parts, then the 
steady-state solution to the Riccati equation is given directly as 

lim P(t) = W22 W-12"^ 

tl-KO 


The program CNTNREG evaluates the solution equations from an input value t-| 
to time zero or until P(t) converges to a steady-state value, whichever 
occurs first. Numerically, convergence is assumed to have occurred when the 
improvement in the element of largest magnitude (measured relatively if the 
magnitude is less than unity, and absolutely otherwise) is less than the 
parameter RICTCV found in the COMMON block CONV of subroutine RDTITL. CNTNREG 
does not evaluate the optimal value of the criterion function or print a 
response trajectory. Options are provided to compute directly the steady- 
state solution and to compute Z, W, S, A, and only without comput- 

ing any of the P(t) values. 

Source of softwar e: Ernest S. Armstrong, LaRC 

Calling sequence: CALL CNTNREG(A, NA,B,NB,H,NH,Q,NQ,R, NR, Z,W, LAMBDA, S,F,NF, 

P , NP , T , lOP , IDENT , DUMMY ) 


Input arguments: 


A, B, H, 

Q, R 

NA, NB, NH, 
NQ, NR 


P 


NP 


Matrices packed by columns in one-dimensional arrays. 

Upon normal return, all except Q are not destroyed. 

Two-dimensional vectors giving the number of rows and col- 
umns in the respective matrices; for example, 

NA ( 1 ) = Number of rows of A 
NA(2) = Number of columns of A 

Not destroyed upon normal return except for NQ replaced 
by NA 

Symmetric nonnegative definite matrix packed in a one- 
dimensional array. On input, P contains the matrix P-|. 
Afterward, intermediate values of P(t) are stored in the 
array P . 

Two-dimensional vector giving the number of rows and columns 
in P: 

NP ( 1 ) = Number of rows of P 
NP(2) = Number of columns of P 
Not destroyed upon normal return 
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T Two-dimensional vector: 

T(1) = Final time t-| 

T(2) = Increment of time for transient solution 
computation 

The final time t-| must be expressed as an integer mul- 
tiple of T(2). The vector T is not required if only 
the steady-state solution of P(t) is required but must 
still appear as an argument of the calling sequence; not 
destroyed upon normal return. Set T(1) negative if 
only Z, W, S, A, and required without 

any P(t) computation. 

lOP Three-dimensional option vector: 

I0P(1) =0 Do not print results, nor compute W”^ZW. 

Otherwise Print input, Z, X^, W, A, W^ZW, Wii, 

Wl2> W21, W22, e-AT(2) and 

values of F and P at time zero or 
steady state, whichever comes first. 

I0P(2) =0 Do not print at intermediate times (multi- 
ples of T(2)) between T(1) and zero or 
steady state. 

Otherwise Regardless of the printing specified by 
I0P(1), print values of P and F at 
intermediate times. 

The parameter I0P(2) is not required if only a steady- 
state solution is required. 

I0P(3) = 0 Compute transient solutions P(t) and 
F(t). 

Otherwise Compute only steady-state values of P 
and F. 

IDENT Logical variable: 

TRUE If H is an identity matrix 

FALSE Otherwise 

For IDENT = TRUE, no data need be input in H but the 
related arguments should still appear in the calling 
sequence . 

DUMMY Vector of working space for computations dimensioned at 

least: 

8n^ + l8n + nr for I0P(3) ^ 0 

9n^ + 17n + nr for I0P(3) = 0 

Output arguments : 

Q Upon normal return, Q is replaced by H'QH- The array Q 

must be dimensioned at least n^. 
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Z, W, LAMBDA, S Matrices packed by columns in one-dimensional arrays dimen- 
sioned at least 4n^, n^, and n^, respectively. 

Upon normal return, these arrays contain their theoretical 
counterparts. The array LAMBDA (A) is declared real. 

F Matrix packed by columns in a one-dimensional array of 

dimension at least nr. Upon normal return, F contains 
the value of F(t) at time zero or the time at which 
steady-state convergence numerically occurs. 

NF Two-dimensional matrix giving, upon normal return, the 

number of rows and columns of F: 

NF(1) = NB(2) 

NF(2) = NA(1) 

P Upon normal return, P contains the value of P(t) at 

time zero or the time at which steady-state convergence 
numerically occurs. 

DUMMY Upon normal return, the first nr elements contain the 

matrix R~^B' , the next 4n^ contain the submatrices 
wn, W 21 , W 12 , and W 22 stored in block column form 

and, if a transient solution is computed, the next n^ 
contain the matrix 


COMMON block ; CONV 
Error messages ; 

(1) If the computation of R“^B’ fails, the message "IN CNTNREG, THE 

SUBROUTINE SYMPDS HAS FOUND THE MATRIX R NOT SYMMETRIC POSITIVE 
DEFINITE" is printed, and the program is returned to the calling 
point. 

(2) If the eigenvalue/eigenvector computation for Z fails, either the 

message "IN CNTNREG, THE EIGENVALUE OF Z HAS NOT BEEN FOUND 

AFTER 30 ITERATIONS" or "IN CNTNREG, EIGEN FAILED TO COMPUTE THE 

EIGENVECTOR OF Z" is printed, and the program is returned to the 
calling point. 

(3) If the computation for W”^ZW fails, the message "IN CNTNREG, GELIM HAS 

FOUND THE REORDERED MATRIX W TO BE SINGULAR" is printed, and compu- 
tation continues. 

(4) If the computation of S for the transient solution fails, the message 

"IN CNTNREG, GELIM HAS FOUND THE MATRIX W21 - P1X W11 TO BE 
SINGULAR" is printed, and the program is returned to the calling 
point. 

(5) If the computation of Wi2~^ fails in the steady-state case, the message 

"IN CNTNREG, GELIM HAS FOUND THE MATRIX W12 TO BE SINGULAR" is 
printed, and the program is returned to the calling point. 
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(6) If at any time, computation of the matrix inverse in the P(t) computa~ 

tion fails, the message "IN CNTNREG AT TIME P CANNOT BE COMPUTED 

DUE TO MATRIX SINGULARITY IN GELIM" is printed, and the program is 
returned to the calling point. 

Field length ; 3046 octal words (1574 decimal) 

Subroutines employed by CNTNREG: EQUATE, TRANP, SYMPDS, SCALE, MULT, JUXTR, 

EIGEN, GELIM, PRNT, SUBT, EXPSER, MAXEL, LNCNT, NULL 

Subroutine employing CNTNREG : ASYMREG 

Comments : Of course, if w(t) = 0 for all t, then the algorithm of CNTNREG 
produces the solution to the deterministic optimal linear output regulator 
problem. In this case, the system response x(t) to the steady-state 
optimal control law 


u(t) = -R-'>B'W 22 Wi 2 "’' x(t) 


is given by 


x(t) = W-| 2 e“'/^^W-] 2 ”^xo 
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Subroutine RICTNWT 


Description ; The purpose of RICTNWT is to solve either the continuous or dis- 
crete steady-state Riccati equation by the Newton algorithms described by 
Kleinman (ref- 9) and Hewer (ref. 10). 

For the continuous case, the algebraic Riccati equation is 


PA + A»P + H*QH - PBR-1b'P = 0 


where the constant matrices A, B, H, Q = Q' ^ 0, and R = R* ^ 0 are of 
dimension n x n, n x r (r ^ n), m x n (m ^ n), m x m, and r x r, respec- 
tively. Applying Newton's equation-solving algorithm (ref. 3) to solve the 
continuous P equation leads to the sequence 


0 = 4>k'Pk + Pk<l>k + H'QH + Fk'RFk 

4>k = A - BFk ) (k = 1,2,...) 

Fk = B-'B'Pk.i 


Kleinman (ref. 9) and Sandell (ref. 2?) established that if the (A,B) pair 
is stabilizable and the (A,D) pair (with D such that D*D = H'QH) is 
detectable (ref. 4), the sequence Pk = Pk' ^ ® (k = 0,1,...) converges to 
the correct P = P' ^ 0 when F-| is chosen so that (A - BF-|) is asymp- 
totically stable in the continuous sense. In RICTNWT, the option is provided 
to solve the continuous Liapunov equation for Pk by either of the BILIN or 
BARSTW subroutines. 

For the discrete case, the steady-state Riccati equation can be written as 


P = 4>*P4) + F'RF + H'QH 


with 

()) = A - BF 

F B (R + B'PB)-''B’PA 

and A, B, H, Q, and R as previously defined. Applying Newton's 
algorithm gives the sequence 
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Pk = 4>k'Pk<t>k + Fk'RFk + h*qh 


Fk = (R + B'Pk_iB)-‘>B'Pk_iA > 


(k H 1,2,...) 


4>k = A - BFk 


It can be established that if the (A,B) pair is stabilizable and the 
(A,D) pair (with D such that D*D = H'QH) is detectable, the sequence 
Pk = Pk' ^ = 0,1,...) converges to the correct P = P* S 0 when Fi 

is chosen so that (A - BF-j) is asymptotically stable in the discrete sense. 
In RICTNWT, the discrete Liapunov equation in Pk is solved by subroutine 
SUM. 


In both the continuous and discrete cases, RICTNWT assumes that an input Fi 
is provided such that (A - F-|) is asymptotically stable in the appropriate 
sense. If A is already stable, F-| = 0 suffices. Otherwise, the subrou- 
tines CSTAB and DSTAB are available for computing initializing Fi matrices. 
RICTNWT repetitively evaluates the solution sequence equations for up to 
100 iterations or until P^ converges to P. Numerically, convergence is 
assumed to have occurred when the improvement in the element of largest magni- 
tude (measured relatively if the magnitude is less than unity, and absolutely 
otherwise) is less than the parameter RICTCV found in COMMON block CONV of 
subroutine RDTITL. Parameters for use in subroutine BARSTW are set internally 
through the COMMON block TOL of RDTITL. 

Source of software ; Ernest S. Armstrong, LaRC 

Calling sequence: CALL RICTNWT(A,NA,B,NB,H,NH,Q,NQ,R,NR,F,NF,P,NP,IOP,IDENT, 

DISC, FNULL, DUMMY) 

Input arguments: 


A, B, H, Matrices packed by columns into one-dimensional arrays. Upon 

Q, R normal return, all except Q are not destroyed. 


NA, NB, NH, Two-dimensional vectors giving the number of rows and columns 
NQ, NR of the respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 

Not destroyed upon normal return except for NQ replaced 
by NA 


F Matrix packed by columns into a one-dimensional array of size 

at least nr. If the matrix A is not asymptotically 
stable, the input F causes (A - BF) to be asymptotically 
stable. For A asymptotically stable, no data for F are 
required. 
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NF 


lOP 


Two-dimensional vector giving the number of rows and columns 
of F; 

NF(1) = r 
NF(2) = n 

Not required as input if data for F are not input 


Three-dimensional option vector: 

I0P(1) =0 Do not print results. 

Otherwise Print input data, H'QH, and final values 
of F and P. 


I0P(2) =0 Do not print at each iteration. 

Otherwise Regardless of printing specified by I0P(1), 
print iteration count and value of P. 


I0P(3) = 0 Use subroutine BARSTW to solve the continuous 
Liapunov equation. 

Otherwise Use subroutine BILIN. 

I0P(3) is not required if the discrete Riccati equation is 
to be solved. 


IDENT Logical variable : 

TRUE If H is an identity matrix 

FALSE Otherwise 

If H is an identity matrix, no data need be input for H, 
but H and NH must still appear as arguments of the 
calling sequence. 

DISC Logical variable : 

TRUE If the digital version is solved 

FALSE For the continuous version 


FNULL Logical variable: 

TRUE If F = 0 

FALSE Otherwise 


DUMMY 

Vector of working space 
5n2 

for 

if 

computations 
DISC = TRUE 

dimensioned at 

least : 


5n^ + n(r + 4) + 2 

if 

DISC = FALSE 

and I0P(3) = 

0 


Tn^ + n(r + 2) 

if 

DISC = FALSE 

and I0P(3) t 

0 


Output arguments : 

Q Upon normal return, Q is replaced by H'QH. Q must be 

dimensioned at least n^. 


F 


The value of 
occurs, F 


Fjj at the stage k of return. If convergence 
contains the steady-state gain matrix. 
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P Matrix packed by columns in a one-dimensional array of 

dimension at least n^; the value of at the stage k 

of return. If convergence occurs, P contains the steady- 
state Riccati equation solution. 

NP, NF Two-dimensional vectors giving the number of rows and columns 

in P and F: upon normal return, 

NP = NA 
NF(1) = r 
NF(2) = n 

COMMON blocks ; TOL, CONV 

Error messages : 

(1) If the computation of either R“^ or (R + B'PjjB)""' (k = 0,1,...) 

fails, the message "IN RICTNWT, A MATRIX WHICH IS NOT SYMMETRIC 
POSITIVE DEFINITE HAS BEEN SUBMITTED TO SYMPDS" is printed, and the 
program is returned to the calling point. 

(2) If the iteration count exceeds 100, the message "THE SUBROUTINE RICTNWT 

HAS EXCEEDED 100 ITERATIONS WITHOUT CONVERGENCE" is printed, I0P(1) is 
set to 1 , P and F are printed, and the program is returned to the 
calling point. 

Field length ; 2415 octal words (1293 decimal) 

Subroutines employed by RICTNWT: LNCNT, PRNT, MULT, TRANP, EQUATE, SYMPDS, 

scale', subt, add, barstw, bilin, maxel, sum 

Subrout ine employing RICTNWT : ASYMREG 

Comments: None 
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Subroutine ASYMREG 


Description ; The purpose of ASYMREG is to solve either the continuous or dis- 
Crete time-invariant asymptotic linear optimal output regulator problem with 
noise-free measurements. For the continuous-time deterministic case, the 
control function u(t) is chosen to minimize the criterion function, 


J = 


lim 

t-i->” 



[y'(t) Q y(t) + u'(t) R u(t^ dt 


subject to 

x(t) = A x( t) + B u(t) 


y(t) = H x(t) 


where x(0) = xq is given and A, B, H, Q=Q'^0, and R = R' > 0 are 
constant matrices of dimension n x n, n x r (r S n), m x n (m ^ n), 
m X m, and r x r, respectively. If the (A,B) pair is stabilizable and the 
(A,D) pair (with D'D = H'QH) is detectable, a solution to the optimal con- 
trol problem exists and is given by 


u(t) = -F x(t) 


where 

' F = R-^B'? 

PA + A'P + H'QH - PBR-Ib'P = 0 


with the criterion function taking the value 


x'(0) P x(0) 


The same control law satisfies the stochastic continuous version of the 
problem (ref. 4) in which the criterion function is 


J = lim E[y'(t) Q y(t) + u'(t) R u(t)3 

t->-co 
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where E is expected value, and the state equation is 


x(t) = A x(t) + B u(t) + w(t) 


where w(t) is a zero-mean Gaussian white noise with intensity V. The 
optimal value of the stochastic criterion function is 


trace (PV) 


For the deterministic discrete-time case, the control function u(i) 
(i = 0,1,...) is chosen to minimize the criterion function 


J 


lim 

N^oo 


N-1 



[y*(i+D Q y(i+1) + u'(i) 


i=0 



subject to 


x(i+1) = A x(i) + B u(i) 


y(i) = H x(i) 


where x(0) = Xq is given and A, B, H, Q, and R are as previously 
defined. If the (A,B) pair is stabilizable and the (A,D) pair (with 
D'D = H’QH) is detectable, a solution to the optimal control exists and is 
given by 


u(i) = -F x(i) 


where 

F = (R + B»PB)-''B'PA 
P = (j>'P(|) + F'RF + H’QH 
(j) = A - BF 


89 



with the criterion function taking the value 


x*(0) P x(0) 


The same control law satisfies the stochastic discrete-time version of the 
problem (ref. 4) in which the criterion function is 


1 

J = lim - E 
lf*CO N 


N-1 

^ [y’(i+D Q y(i+1) 
i=0 


+ u’(i) 



and the state equation is 

x(i+1) = A x(i) + B u(i) + w(i) 


where w(i) (i = 0,1,...) is a sequence of zero-mean stochastic variables 
with variance matrix V. The optimal value of the stochastic criterion func- 
tion is 


trace (PV) 


ASYMREG does not evaluate the optimal values of the performance criteria. 
Therefore, no V data are input. The option of solving the appropriate 
steady-state Riccati equation using either of the subroutines DISCREG, 
CNTNREG, or RICTNWT is provided. In addition, the residual error in the 
Riccati equation, the eigenvalues of P, the closed-loop response matrix 
(A - BF) , and the eigenvalues of (A - BF) can be computed. If RICTNWT is 
selected and if the matrix A is not asymptotically stable or it is unknown 
whether A is asymptotically stable, the option is provided to test the 
relative stability of A and compute, if necessary, a stabilizing gain to 
initialize the Newton process . 

Source of software ; Ernest S. Armstrong, LaRC 

Calling sequence ; CALL ASYMREG(A,NA,B,NB,H,NH,Q,NQ,R,NR,F,NF,P,NP,IDENT,DISC, 
NEWT , STABLE , FNULL , ALPHA , lOP , DUMMY ) 

Input arguments : 

A, B, H, Matrices packed by columns into one-dimensional arrays. Upon 

Q, R normal return, all except Q are not destroyed. 
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NA, KB, 
NQ, NR 


F, NF, 
P, NP 


IDENT 


DISC 


NEWT 


STABLE 


FNULL 


NH, Two-dimensional vectors giving the number of rows and eoliamns 
of the respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 

Not destroyed upon normal return except for NQ replaced 
by NA 

F and P are matrices packed by columns into one-dimensional 
arrays of dimension at least nr and n^, respectively. 

NF and NP are the corresponding two-dimensional vectors 
giving the number of rows and columns of P and F 
Any input requirements depend on and are specified by whether 
DISCREG, CNTNREG, or RICTNWT is employed. 

Logical variable: 

TRUE If H is an identity matrix 

FALSE Otherwise 

If H is an identity matrix, no data need be input for H, 
but H and NH must still appear as arguments of the 
calling sequence. 

Logical variables: 

TRUE If the digital version is solved 
FALSE For the continuous version 

Logical variable: 

TRUE Newton's method (RICTNWT) is used to solve the 

appropriate Riccati equation. 

FALSE Either CNTNREG or DISCREG is used depending upon 
the value of DISC. 

Logical variable. Data for STABLE are not required if 
NEWT = FALSE, but STABLE must still appear as an argument 
of the calling sequence. When NEWT = TRUE: 

STABLE = TRUE If it is known that the matrix (A - BF) 

computed from input data is stable 
relative to an input parameter ALPHA. 
STABLE = FALSE The matrix (A - BF) is evaluated and 

tested for stability relative to ALPHA 
using subroutine TESTSTA. If a stabi- 
lizing gain is required, it is computed 
from subroutine DSTAB or CSTAB. 

Logical variable; data for FNULL are not required if 

NEWT = FALSE, but FNULL must still appear as an argument of 
the calling sequence. When NEWT = TRUE: 

FNULL = TRUE If the input F is a null matrix 
FNULL = FALSE Otherwise 
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alpha Scalar variable. ALPHA is not required if NEWT = FALSE or 

STABLE = TRUE. Otherwise, ALPHA is used in the asymptotic 
stability test of (A - BF) from input data using subrou- 
tine TESTSTA. 


lOP 


Five-dimensional option vector: 

I0P(1), I0P(2), and I0P(3) are the first three elements of 
the lOP vector in DISCREG, CNTNREG, or RICTNWT. If CNTNREG 
is selected, I0P(3) is set internally to a nonzero value. 


I0P(4) =0 Do not compute the Riccati equation residual. 
Otherwise Compute the residual and print. 


I0P(5) = 0 Compute but do not print the eigenvalues 

of P, the matrix (A - BF) , and the eigen^ 
values of (A - BF). 

Otherwise Print these data after computation. 


DUMMY Vector of working space for computations dimensioned at least: 


DISC 

TRUE 

NEWT 

TRUE 

6n^ + 4n + 2 

(STABLE = FALSE) 

TRUE 

FALSE 

5n^ 

4n2 

(STABLE = TRUE) 

FALSE 

FALSE 

17n^ + n(r + 18) 


FALSE 

TRUE 

8n^ + n(r + 1 ) 

(BILIN) 



5n^ +n(r+4) +2 

(BARSTW) 


Output arguments : 


Q 


Upon normal return, Q is replaced by H*QH. The dimension 
of Q must be at least n^. 


F Upon normal return, F contains the steady-state gain matrix 

for either the continuous or discrete problem. 

NF Two-dimensional vector giving, upon normal return, the number 

of rows and columns of F: 

NF(1) = r 
NF(2) = n 


P Upon normal return, P contains the steady-state solution for 

either the continuous or discrete Riccati equation. 

NP Two-dimensional vector giving the number of rows and columns 

of P. Upon normal return, NP = NA. 

STABLE Upon normal return, STABLE = TRUE. 
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DUMMY 


Upon normal return, the first n elements of DUMMY contain 
the eigenvalues of P, the next n^ contain the matrix 
(A - BF) for steady-state F, and the next 2n contain 
the eigenvalues of (A - BF) stored as an n x 2 matrix 
with the real parts as first column and imaginary parts as 
the second. All matrices are packed by columns into one- 
dimensional arrays. 

COMMON blocks ; None 

Error messages ; 

(1) If the stabilizing gain computation for the continuous system fails, the 

message "IN ASYMREG, CSTAB HAS FAILED TO FIND A STABILIZING GAIN 

MATRIX (F) RELATIVE TO / ALPHA = " is printed, and the program 

is returned to the calling point. 

(2) If the stabilizing gain computation for the discrete system fails, the 

message "IN ASYMREG, DSTAB HAS FAILED TO FIND A STABILIZING GAIN 

MATRIX (F) RELATIVE TO / ALPHA = " is printed, and the program 

is returned to the calling point. 

(3) If the eigenvalue computation for P fails, the message "IN ASYMREG, 

THE EIGENVALUE OF P HAS NOT BEEN COMPUTED AFTER 30 ITERATIONS" 

is printed, and the program continues with the computation of (A - BF) 
and its eigenvalues . 

(4) If the eigenvalue computation for (A - BF) fails, the message "IN 

ASYMREG, THE EIGENVALUE OF A-BF HAS NOT BEEN COMPUTED AFTER 

30 ITERATIONS" is printed, and a return is made to the calling point 
if no printing is required. Otherwise, the available information is 
printed before return. 

Field length : 1756 octal words (1006 decimal) 

Subroutines employed by ASYMREG: MULT, SUBT, TESTSTA, CSTAB, SCALE, DISCREG, 

RICTNWT, TRANP, ADD, EQUATE, EIGEN, JUXTC, LNCNT, PRNT, CNTNREG, DSTAB 

Subroutines employing ASYMREG : ASYMFIL, EXPMDFL, IMPMDFL 

Comments ; When using DISCREG to solve the steady-state discrete Riccati equa- 
tion, the entry I0P(3) = N must be set sufficiently large for steady-state 
convergence to occur. 
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Subroutine ASYMFIL 


Description ; The purpose of ASYMFIL is to solve either the continuous or dis- 
crete time-invariant asymptotic optimal Kalman-Bucy filter problem (ref. 4). 

For the continuous case, the state equation is given by 
x(t) = A x(t) + G n(t) 


with output 

y(t) = H x(t) + m(t) 


where A, G, and H are constant matrices of dimension n x n, n x m 
(m S n) , and r x n (r ^ n) , respectively. The process noise n(t) is a 
zero-mean Gaussian white-noise process with intensity for the covariance 
given by Q = Q' ^ 0. The measurement noise m(t) is also a zero-mean 
Gaussian white-noise process with intensity matrix R = R' > 0. The pro- 
cesses n(t) and m(t) along with the Gaussian x(to) are mutually uncor- 
related. The optimal filter problem is to construct an estimate x(t) of 
x(t) operating over |^to,tJ such that the quantity, 

J = lim E[e'(t) W e(t)] 
to^-oo 

is minimized where E denotes expected value and 
e(t) = x(t) - x(t) 


and the constant n x n matrix satisfies 


W = W ^ 0 


When the pair (A',H') is stabilizable and the (A*,D') pair (where 
DD' = GQG') is detectable, the solution to the asymptotic optimal observer 
problem exists and x(t) is given by 

x(t) = A x(t) + F[y(t) - H x(t)] 
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»rtiere 


x(to) = E[x(to)] 


with filter gain 


F = PH'R-I 


and P satisfying 


0 = AP + PA' + GQG' - PH'R-IhP 


The matrix P = P' ^ 0 represents the (constant) steady-state variance 
matrix of the reconstruction error e(t); that is, 

lim E[e(t) e'(t)] = P 
to-^-“ 


Also, 


J = lim E[e'(t) W e(t)^ = trace (PW) 

tQ-^-CO 


For the discrete case. 


x(i+1) = A x(i) + G n(i) 


with output 


y(i) = H x(i) + m(i) 


where A, G, and H are as previously defined. The process noise n(i) 

(i = io ,io+1 > • • • ) is a sequence of zero-mean Gaussian white-noise stochastic 
variables with variance matrix Q = Q' S 0. The measurement noise m(i) 

(i = io,io+1»**») is also a zero-mean Gaussian white-noise (ref. 4) sequence 
with variance R = R' > 0. The processes n(i) and m(i) along with the 
Gaussian x(ig) are mutually uncorrelated. The optimal estimator problem 
considered here (technically a prediction problem) is to construct an esti- 
mate x(i) of x(i) from knowledge of y(io) ,y(io+1 ) » • • • »y(i-1 ) such that 
the quantity. 
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J = lim E[e*(i) W e(i^ 

iQ^-OO 

is minimized, tdiere 

e(i) = x(i) - x(i) 


with W as previously defined. When the (A',H') pair is stabilizable and 
the (A',D’) pair (with DD* = GQG') is detectable, the solution to the 
asymptotic optimal observer problem exists and x(i) is given by 

x(i+1) = A x(i) + pjyCi) - H x(i)J 


where 


x(io) = E[x(io)J 


with filter gain 


F = APH'(R + 


and P satisfying 


P = W + FRF’ + GQG’ 


for 


4i = A - FH 


The matrix P = P' ^ 0 represents the (constant) steady-state variance 
matrix of the reconstruction error e(i); that is. 


lim E[e(i) e’(i)] = P 


Also, 


J = lim E[e'(i) W e(i)] = trace (PW) 

ig-^-CO 
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Computation for both the discrete and continuous versions of the foregoing 
optimal filter problems is performed using duality theory (ref. M) and the 
regulator subroutine ASYMREG. Thus, the user has the option of solving the 
appropriate steady-state covariance equations in P by either of the sub- 
routines DISCREG, CNTNREG, or RICTNWT. No computations are performed which 
involve the matrix W; hence, no data for W are required. 

Sou rce of software ; Ernest S. Armstrong, LaRC 

Calling sequence ; CALL ASYMFIL(A,NA,G,NG,H,NH,Q,NQ,R,NR,F,NF,P,NP,IDENT,DISC, 
NEWT , STABLE , FNULL .ALPHA , lOP , DUMMY ) 

Input arguments ; 

A, G, H, Matrices packed by columns in one-dimensional arrays. Upon 

Q, R normal return, all except Q are not destroyed. 

NA, NG, NH, Two-dimensional vectors giving the number of rows and columns 
NQ, NR of respective matrices; for example, 

NA(1) = Number of rows of A 
NA(2) = Number of columns of A 

Not destroyed upon normal return except for NQ replaced 
by NA 

F, NF, F and P are matrices packed by columns into one-dimensional 

P, NP arrays of dimension at least rn and n^, respectively. 

NF and NP are the corresponding two-dimensional vectors 
giving the number of rows and columns in F and P. 

Any input requirements depend on and are specified by whether 
DISCREG, CNTNREG, or RICTNWT is employed by ASYMREG. 

When F data are input, it should be kept in mind that, 
because of the use of duality theory, ASYMREG will treat A* 
as A and H' as B. An initial stabilizing gain F for 
use in RICTNWT must cause (A' - H'F) to be asymptotically 
stable and be appropriately dimensioned. 

IDENT Logical variable: 

TRUE If G is an identity matrix 
FALSE Otherwise 

If G is an identity matrix, no data need be input for G, 

but G and NG must still appear as arguments of the call- 
ing sequence. 

DISC Logical variable: 

TRUE If the digital version is solved 
FALSE For the continuous version 

NEWT Logical variable: 

TRUE RICTNWT is used to solve the appropriate steady- 

state covaricihce equation 

FALSE Either CNTNREG or DISCREG is used depending upon 
the value of DISC. 
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STABLE 


FNULL 


ALPHA 


lOP 


Logical variable. Data for STABLE are not required if 

NEWT = FALSE, but STABLE must still appear as an argiament 
of the calling sequence. When NEWT = TRUE: 

STABLE = TRUE If it is known that the matrix (A' - H'F) 

computed from input data is stable 
relative to an input parameter ALPHA. 
STABLE = FALSE (A' - H'F) is computed and tested 

(within ASYMREG) for stability relative 
to ALPHA, 


Logical variable. Data for FNULL are not required if 

NEWT = FALSE, but FNULL must still appear as an argument 
of the calling sequence. When NEWT = TRUE: 

FNULL = TRUE If the input F is a null matrix 
FNULL = FALSE Otherwise 


SCALAR VARIABLE. ALPHA is not required if NEWT = FALSE or 
STABLE = TRUE. Otherwise, ALPHA is used in the asymptotic 
stability test of (A' - H’F) from input data. 


Five-dimensional option vector: 

I0P(1) =0 Do not print results. 

Otherwise Print input, GQG', F, P, the eigenvalues 

of P, the matrix (A - FH) , and the eigen- 
values of (A - FH) . 


I0P(2), I0P(3), I0P(4), and I0P(5) are the first four 
elements of the lOP vector of ASYMREG. 


DUMMY Vector of working space for computations, dimensioned at 


least: 




DISC 

NEWT 



TRUE 

TRUE 

6n^ + 4n + 2 

(STABLE = FALSE) 



5n2 

(STABLE = TRUE) 

TRUE 

FALSE 

4n^ 


FALSE 

FALSE 

17n^ + n(r + 18) 


FALSE 

TRUE 

8n^ + n(r +1) 

(BILIN) 



5n^ +n(r+M) +2 

(BARSTW) 


Output arguments ; 

Q Upon normal return, Q is replaced by GQG'. The dimension 

of Q must be at least n^. 

F Upon normal return, F contains the filter gain for either 

the continuous or the discrete problem. 

NF Two-dimensional vector giving, upon normal return, the 

number of rows and columns of F: 

NF(1) = n 
NF(2) = r 
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Upon normal return, P contains the steady-state covariance 
matrix. 

Two-dimensional vector giving the number of rows and columns 
of P. Upon normal return, NP = NA. 

Upon normal return, the first n elements of DUMMY contain 
the eigenvalues of P, the next n^ contain the matrix 
(A - FH) for filter gain F, and the next 2n contain 
the eigenvalues of (A - FH) stored as an n x 2 matrix 
with the real parts as first column and imaginary parts as 
the second. All matrices are packed by columns into one- 
dimensional arrays. 

COMMON blocks ; None 

Error messages ; None directly from ASYMFIL 
Field length ; 713 octal words (459 decimal) 

Sub routines employed by ASYMFIL ; LNCNT, PENT, TRANP, EQUATE, ASYMREG 
Subrouti nes employing ASYMFIL : None 

Comments ; Cases in which G, with dimension n x m, has m ^ n can be treated 
as follows. Compute GQG* externally and, using subroutine FACTOR, find a 
q X n matrix D (q < n) such that 


P 

NP 

DUMMY 


GQG' = D'D 


Then apply ASYMFIL with G replaced by D' and Q = Ij^. 

Extensions of the basic optimal filter problem considered in ASYMFIL (such 
as colored noise, singular R matrices, and correlated process and measure- 
ment noise) can be found in the literature (e.g., ref. 4). Generally, each 
extension can be solved by an appropriate combination of ASYMFIL with other 
subroutines of the ORACLS program. The transient Kalman-Bucy filter problem 
in which tQ and ig are finite can be solved by using duality theory and 
the transient regulator solution capability of subroutines CNTNREG and 
DISCREG. 
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Subroutine EXPMDFL 


Description ; The purpose of EXPMDFL is to solve either the continuous or dis- 
crete time- invariant asymptotic explicit (model-in-the-system) model-following 
problem (ref. 28). 

For the continuous case, the state and output equations are given as 


x(t) = A x(t) + B u(t) 


y(t) = H x(t) 


where x(0) = xg is given and the constant matrices A, B, and H are of 
dimension n x n, n x r (r ^ n), and m x n (m ^ n), respectively. The 
control function u(t) is required to minimize 

J = lim I^J* ^ [e'(t) Q e(t) + u(t») R u(t)j dtj 

where 

e(t) = y(t) - ym(t) 

Xm ^ ^ ^ ” Hjn ^m ^ ^ ^ 

and 


*m^^) “ ^m ^m^^^ 


where Xm(0) = x^ is given. The constant matrices and A^ have 

dimension m x 8, (m ^ S,) and i x Z, respectively. Also, Q = Q' S 0 
and R = R* > 0. The optimization of the performance index causes the 
output y(t) of the state to track the output yoj(t) of a prescribed model. 
After substituting e(t) into the performance index, the model-following 
problem can be transformed into choosing u(t) to minimize 


J = lim j [x'(t) Q x(t) + u'(t) R u(t)^ dt 

ti^ 

with 

x(t) = A x(t) + B u(t) 
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where 


A = 


B 


B = 




Q = 


H’QH 

-H'QH, 

-Hm'QH 

Hffl’QH, 


and 


X 


X 


X 


m 


This transformed problem_^can be solved directly using optimal ^linear regu- 
lator theory.^ If the (a,b) pair is stabilizable and the (a,d) pair 
(with D'D = q) is detectable, the solution exists and is given by 


u(t) = -F x(t) = -F-ii x(t) - F -|2 Xn,(t) 


Computationally, it is inefficient to work with the composite LA,BJ system 
directly. If the steady-state Riccati equation is formed and A, B, and 
Q are substituted, it readily follows (ref. 28) that 


Fii = R-Ib'Pii 

with Pii = Pii' = 0 satisfying 

PllA + A»Pii - PiiBR-IB'Pii + H'QH = 0 
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and 

F-|2 = R"''b*Pi2 


with P-]2 satisfying 

P-|2An, + (A - BFii)’Pi2 = H-QH^ 

The computation of (Fii»Fl2^ thus separates into two parts: 

(1) Evaluate the feedback gain F-| -) on the state x by solving a reduced- 
order optimal regulator problem of the form 


x(t) = A x(t) + B v(t) 


y(t) = H x(t) 


min 

v(t) 



[y'(t) Q y(t) + v'(t) R v*(t)] dt 


leading to 


v(t) = -F-|i x(t) 

(2) Using F-]-] from step (1), compute the feedforward gain F-|2 the 
model Xjn from the linear equations 

Pl2Am + (A - BFii)'Pi2 = H'QH„ 

RF 12 = ® 2 

For the discrete case, the state and output equations are given as 
x(i+1) = A x(i) + B u(i) 
y(i) = H x(i) 
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with A, B, and H as previously defined. The control sequence u(i) 
(i = 0,1,..., N-1) is required to minimize 


where 


J 



N-1 
i=0 


e*(i+1) Q e(i+1) + u'(i) 



e(i) = y(i) - ymCi) 

ym(i) “ Njjj X|j](i) 
Xm(i+1) = Am Xn,(i) 


with Q, R, Hm, and Am as previously defined. As in the continuous^case, 
the discrete model-following problem can be solved in terms of a (a,B,Q,r) 
optimal regulator formulation, but a simplified computational algorithm also 
exists (ref. 15): 

(1) Compute a feedback gain F-| •) on the state x by solving the reduced- 
order optimal regulator problem 


x(i+1) = A x(i) + B v(i) 


y(i) = H x(i) 


N-1 


in / lim \ []y'(i+1) Qy(i+1) +v’(i) R v(i)J 

( 1 ) 1 1^0 J 


min 

V 


leading to 


v(i) = -Fii x(i) 


(2) Using (P-|i,Fii) from step (1), compute a feedforward gain F -\2 on the 
model state Xm from the linear equations, 


Pl2 = (A - BFii)'Pi2Am - H*QHm 


(B’PiiB + R)Fi 2 = B’Pi2Am 
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The complete optimal model-following control law is then given by 
u(i) = -Fii x(i) - F^2 Xin(i) 


EXPMDFL solves both the continuous and discrete versions of the explicit 
model-following problem through the simplified computational approach of 
solving a reduced-order regulator program for F-| i and a set of linear equa- 
tions for Fi 2- Subroutine ASYMREG is used to solve the reduced-order regu- 
lator problem, thus giving the user the choice of solving the appropriate 
steady-state Riccati equation for P-|i by either of the subroutines DISCREG, 
CNTNREG, or RICTNWT. The subroutine BARSTW is used to solve the P-|2 equa- 
tion in the continuous case and subroutine SUM in the discrete case. Final 
F -|2 computation in both cases employs subroutine SYMPDS. Computational 
parameters for BARSTW are set internally through the COMMON block TOL of 
subroutine RDTITL. An option is provided to bypass the (Pii,Fii) computa- 
tion and, using input (P'|i,Fii) data, proceed directly to the computation 
of ‘ 

Source of softw are; Ernest S. Armstrong, LaRC 

Calling seq uence; CALL EXPMDFL(A,NA,B,NB,H,NH,AM,NAM,HM,NHM,Q,NQ,R,NR,F,NF, 

P , NP , HIDENT , HMIDENT , DISC , NEWT , STABLE , FNULL , ALPHA , lOP , DUMMY ) 


Input argumen ts : 

A, B, H, AM, Matrices packed by columns in one-dimensional arrays; not 

HM, Q, R destroyed upon normal return 


NA, NB, NH, Two-dimensional vectors giving the number of rows and 

NAM, NHM, columns of the respective matrices; for example, 

NQ, NR NA(1) = Number of rows of A 

NA(2) = Number of columns of A 
Not destroyed upon normal return 


F 


NF 


Matrix packed by columns into an array dimensioned at 
least r(n + 1 ). The first rn elements of F are 
treated as an input matrix F to ASYMREG for solving 
the reduced-order Riccati equation in Pi •] . If no input 
F data are required by ASYMREG, no data need be input 
for F in EXPMDFL unless the reduced-order Riccati 
computation for (Pn,Fii) is to be bypassed. In this 
case, the first rn elements of F should contain a 
matrix Fn to be used in the computation of (A - BFn). 

Two-dimensional vector giving, if required, the number of 
rows and columns of F used by ASYMREG or Fn; 

NF(1) = r 
NF(2) B n 
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p 


NP 


HIDENT 


HMIDENT 


DISC 


Matrix packed by columns into an array dimensioned at 
least n(n + Z) . The first n^ elements of P are 
treated as an input matrix P to ASYMREG when used to 
solve the reduced-order Riccati equation for Pii- If 
no input P data are required by ASYMREG, no data need 
be input for P in EXPMDFL unless the reduced-order 
Riccati computation for (Pi-|,F-|i) is to be bypassed in 
the discrete case. In the discrete case, the first n^ 
elements must contain a matrix to be used for Pii in 
the computation of the discrete F-i2* 

Two-dimensional vector giving, if required, the number of 
rows and columns of P used by ASYMREG or Pi -| ; 

NP(1) = n 
NP(2) = n 

Logical variable: 

TRUE If H is an identity matrix 

FALSE Otherwise 

If H is an identity matrix, no data need be input 
for H, but H and NH must still appear as arguments 
of the calling sequence. 

Logical variable: 

TRUE If Hm (HM) is an identity matrix 
FALSE Otherwise 

If Hm is an identity matrix, no data need be input 
for HM, but HM and NHM must still appear as argu- 
ments of the calling sequence. 

Logical variable: 

TRUE If the discrete version is solved 

FALSE For the continuous version 


NEWT, STABLE, 
FNULL, ALPHA 


lOP 


Variables whose input values are determined by the choice 
of method to solve the steady-state Riccati equation for 
(Pl 1 ,Fii) called for in ASYMREG; not required if the 
(Pl 1 ,Fii) computation is bypassed but still must appear 
as arguments of the calling sequence 


Five-dimensional option vector: 

I0P(1) =0 Do not print results. 

Otherwise Print input and computed results. 


I0P(2) =0 Do not compute (Pii,Fn) but use 

input P and F as these variables. 
Otherwise Compute Pn and Fn through ASYMREG. 

I0P(3)> I0P(4), and I0P(5) are first three elements of 
the lOP vector in ASYMREG; not required if I0P(2) = 0. 


i 
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DUMMY Vector of working space for computations, dimensioned at 

least: 

I0P(2) DISC 

0 TRUE 

0 FALSE 

Nonzero TRUE or 

FALSE 


Output arguments : 

F Upon normal return, F = [Fii,Fi2j packed by columns in 

a one-dimensional array 

NF Upon normal return, 

NF(1) = r 
NF(2) = n + il 

P For I0P(2) ^0, P = upon normal return. 

If I0P(2) = 0, the first nJl elements contain the 
matrix P-i2* In both cases, P is packed by columns 
into a one-dimensional array. 

NP Upon normal return, 

NP(1) = n 

NP(2) = n + S- (I0P(2) ^ 0) 

or 

NP(2) □ £ (I0P(2) B 0) 

COMMON block : TOL 

Error message: If either R or (B'PnB+R) fails to be positive definite, 

the message "IN EXPMDFL, THE COEFFICIENT MATRIX FOR SYMPDS IS NOT SYMMETRIC 
POSITIVE DEFINITE" is printed, and the program is returned to the calling 
point . 

Field length : 1516 octal words (846 decimal) 

Subroutines employed by EXPMDFL: LNCNT, PRNT, EQUATE, ASYMREG, MULT, SUBT, 

TRANP, BARSTW, SUM, ADD, SYMPDS, SCALE 

Subroutines employing EXPMDFL: None 


3n^ + £.2 + max(n2,£n) 

3n2 + 4n + 2 + max(n2,£,n) 

Either n^ plus the DUMMY require- 
ment of ASYMREG or the preceding 
requirements for I0P(2) = 0, 
whichever is larger 
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Comments ; If the model dynamics in both the continuous and discrete cases 
contain a direct link to the control, then 


B = 



(Bn i 0) 


and the computation will not generally decouple into a simglified algorithm. 
In this case, ASYMREG can be applied to the composite (A,B,Q,R) system 
directly. For the transient case where N and are finite, the time- 

varying explicit model-following solution can be obtained using the composite 
(A,S,Ci,R) system and subroutine DISCREG or CNTNREG. 
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Subroutine IMPMDFL 


Description ; The purpose of IMPMDFL is to solve either the continuous or dis- 
crete time-invariant asymptotic implicit (model-in-the-performcince-index) 
model-following problem. 

For the continuous case, the state and output equations are given as 


x(t) = A x(t) + B u(t) 


y(t) = H x(t) 


where x(0) = xq is given and the constant matrices A, B, and H are of 
dimension n x n, n x r (r = n) , and m x n (m S n) , respectively. The con- 
trol function u(t) is required to minimize 


J = 



Q e(t) + u'(t) R u(t)J dt 


} 


where 

e(t) = y(t) - A„j y(t) - Bj„ u(t) 

The constant matrices Ajn and B^ have dimension m x m and m x r, respec- 
tively. Also, Q = Q’ ^ 0 and R = R* >0. The philosophy of implicit model 
following is to cause the output y(t) of the system state to behave similar 
to a model state with dynamics 


Xm(t) = Ajj] Xjjj(t) + Bjjj u(t) 


o 

where Xj^CO) = Xj„ is given. This is done by forcing, in a weighted least 
squares sense, y(t) to satisfy the model dynamic equation. Substituting 
for e(t) in the performance index reduces the problem to choosing u(t) 
to minimize 


J 


lim 



Q x(t) + x'(t) W u(t) + u'(t) R u(t)J 
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where 


Q = (HA - AmH)'Q(HA - AnjH) 


W H 2(HA - AnjH)*Q(HB - Bn,) 




and 


R = (HB - Bn,)'Q(HB - B„,) + R 


ni'' 


Applying the control transformation 


u(t) H -Fq x(t) + v(t) 


with 


~_1 W’ 
Fo = R - 


further reduces the problem to choosing v(t) to minimize 


= 11m f r Tx'Ct) Q x(t) + v'(t) R v(t)H dt 


where 


/V ~ W 
Q = Q - - Fo 


A = A - BFq 


and 


x(t) = A x(t) + B v(t) 
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If the (A,B) pair is stabilizable and the (A,D) pair (with D*D = Q) is 
detectable, a solution to the asymptotic implicit model-following problem 
exists and is given by 


v(t) = -F-i x(t) 


PA + A*P + Q - PbF^B'P = 0 


and 


u(t) = -F x(t) 


with 


F — Fq + F'^ 


For the discrete case, the state and output equations are given as 


x(i+1) = A x(i) + B u(i) 


y(i) = H x(i) 


where x(0) = xq is given and A, B, and H are as previously defined. 
The control function u(i) (i = 0,1,..., N-1) is required to minimize 


J = lim 

IJ>oo 



[e'(i) Q e(i) + u*(i) 



where 


e(i) = y(i+1) - Am y(i) - 6^ u(i) 
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The matrices Aj^, Bj^, Q, and R are as previously defined. Substituting 
for e(i) in the performance index reduces the problem to choosing u(i) 
to minimize 


J 


rN-1 

lim I V [x'(i) Q x(l) + x'(i) W u(i) + u’(i) 
"■*" 1 1=0 



Applying the control transformation, 
u(i) = -Fq x(i) + v(i) 


further reduces the problem to choosing v(i) to minimize 

fv 

J = x’(0) Q x(0) + lim < y [x’(i+1) Q x(i+1) + v'(i) R v(i)J 

N-^oo 1 ^ 

L 1=0 > 


with 


x(i+1) = A x(i) + B v(i) 


If the aforementioned stabilizabillty and detectability conditions are 
satisfied 


v(i) = -F-] x(i) 

Fi = (R + B'PB)-''b*P a 
P = <t)'P(t) + Ft' R F-) + Q 


4) = A - BF-i 


and 


u(i) D -F x(i) 
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with 


F = Fq + Fi 


The program IMPMDFL solves the steady-state Riccati equation to generate the 
appropriate P and F-) through use of the subroutine ASYMREG, thus giving 
the user the choice of either subroutine DISCREG, CNTNREG, and RICTNWT to 
solve the final Riccati equation. An option is provided^to bypass the 
(P,F-|) computation and return after computing Q, W, R, Fq, Q, and A. 

Source of software : Ernest S. Armstrong, LaRC 

Calling sequence ; CALL IMPMDFL(A,NA,B,NB,H,NH,AM,NAM,BM,NBM,Q,NQ,R,NR,F,NF, 

P , NP , IDENT , DISC , NEWT , STABLE , FNULL , ALPHA , iOP , DUMMY ) 


Input arguments ; 

A, B, H, AM, Matrices packed by columns in one-dimensional arrays; 

BM, Q, R not destroyed upon normal return 


NA, NB, NH, Two-dimensional vectors giving the number of rows and 

NAM, NBM, columns of respective matrices; for example, 

NQ, NR NA(1) = Number of rows of A 

NA(2) = Number of columns of A 
Not destroyed upon normal return 


F, NF, P, NP F and P are matrices packed by columns into one- 

dimensional arrays of dimension at least rn and n^, 
respectively. 

NF and NP are the corresponding two-dimensional vectors 
giving the number of rows and columns of F and P. 

Any input requirements depend on and are specified by 
whether DISCREG, CNTNREG, or RICTNWT is employed by 
ASYMREG. When F data are input for RICTNWT, it should 
be such that the input (A,B) system is stabilized. 


IDENT Logical variable: 

TRUE If H is an identity matrix 

FALSE Otherwise 

If H is an identity matrix, no data need be input 
for H, but H and NH must still appear as arguments 
of the calling sequence . 

DISC Logical variable: 

TRUE If the discrete version is solved 

FALSE For the continuous version 


NEWT, STABLE, Variables whose input values are determined by the choice 
FNULL, ALPHA of method called for from ASYMREG to solve the steady- 

state Riccati equation for P and F 
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Four-dimensional option vector; 

lOP(l) =0 Do not print results. 

Otherwise Print input and computed results. 

I0P(2), I0P(3), and I0P(4) are the first three elements 
of the lOP vector in ASYMREG. If I0P(2) = -1000, the 
IMPMDFL program returns just prior to employing ASYMREG, 
in which case data for F, NF, P, NP, DISC, NEWT, STABLE, 
FNULL, ALPHA, I0P(3), and I0P(4) need not be entered, 
but these arguments must still appear in the calling 
sequence. 

Vector of working space for computations, dimensioned at 
least 4n^ plus the dimension requirements for DUMMY in 
ASYMREG when (P,F-i) is to be computed. If the (P,F-|) 
computation is bypassed, DUMMY should be dimensioned at 
least 6n2 + r. 


If (P,F-|) is computed, F = Fg + F-j upon normal return. 
For (P,F-|) not computed, the input F is returned. 

All matrices are packed by columns into a one-dimensional 
array. 

Upon normal return, 

NF(1) = r 
NF(2) = n 

If (P,Fi) is computed, P contains, upon normal return, 
the solution to the steady-state Riccati equation defin- 
ing Fi . For (P,F-|) not computed, the input P is 
returned. All matrices are packed by columns into one- 
dimensional arrays. 

Upon normal return, 

NP(1) = n 
NP(2) = n 

when (P,F-|) is computed. Otherwise, the input NP 
is returned. 

Upon normal return, if I0P(2) = -1000, the first n^ ele- 
ments contain the matrix Q and the next nr contain 
the matrix Fg. After 2n^ elements, the next r^ con- 
tain the matrix R. After 3n^, the next r?- contain 
the matrix A. Additionally, for (P,F-|i) computed, 
after 4n^ , the elements of DUMMY in IMPMDFL contain the 
elements of DUMMY returned by ASYMREG. All matrices are 
packed by columns as one-dimensional arrays. 



COMMON blocks: None 


Error messages : None directly from IMPMDFL 

Field length : 1556 octal words (8?8 digital) 

Subroutines employed by IM PMDFL: LNCNT, PRNT, SUBT, EQUATE, PREFIL, ASYMREG, 

ADD, MULT, TRANP, SCALE ' 

Subroutines employing IMPM DFL : None 

Comments : For the case in which t-| or N is finite, set I0P(2) =^-1000^ 

and return with Q, Fq, R, and A in DUMMY. Then, using the (A,B,Q,R) 
system, solve for the transient P( ) and F-|( ) from subroutine DISCREG 
or CNTNREG. The model-following control law is then 


F( ) = Fq + Fi( ) 
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SUPPORTING SUBROUTINES 


Subroutine READ1 

Description ; The purpose of READ1 is to input from cards and print or print 
without card input a single matrix A without employing a header card as 
in subroutine READ. When card input is employed, each row of the matrix 
starts on a new card using format (8F10.2). Printing is performed using 
subroutine PRNT. 

Source of software ; VASP (ref. 16) with modifications by Ernest S. Armstrong, 
LaRC 

Calling sequence ; CALL READ1 (A,NA,NZ,NAM) 

Input arguments ; 

A Matrix packed by columns in a one-dimensional array; required as 

input only if NZ(1) = 0 

NA, NZ Two-dimensional arrays generally giving the number of rows and 
columns of A: 

NA(1) = NZ(1) = Number of rows of A 
NA(2) = NZ(2) = Number of columns of A 

Alternatively, if NZ(1) = 0, the data for (A,NA) are assumed 
to have been previously stored in locations A and NA. Other- 
wise, data for NA are not required as input. 

NAM Hollerith data: a four-character symbol used by PRNT 

Output arguments ; 

A, NA For NZ(1) ^ 0, the matrix A which was read from cards and packed 
by columns into a one-dimensional array; upon normal return, 

NA(1) = NZ(1) 

NA(2) = NZ(2) 

COMMON blocks ; None 

Error message ; If either NZ(1) or NZ(1) x NZ(2) is less than unity, the 

message "ERROR IN READ1 MATRIX HAS NA = , ” is printed, 

and the program is returned to the calling point. 

Field length ; 133 octal words (91 decimal) 

Subroutines employed by READ1 ; LNCNT, PRNT 

Subroutine employing READ1 ; READ 

Comments; None 
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Subroutine BALANC 


Description ; The purpose of BALANC is to balance a real square matrix for the 
calculation of eigenvalues and eigenvectors and isolate the eigenvalues when- 
ever possible. The computational method follows that of Parlett and Reinsch 
found in reference 17. The subroutine BALANC is a translation of the ALGOL 
procedure BALANCE found on pages 315 to 326 of reference 2. 

Source of softw are : LaRC Analysis and Computation Division subprogram library 

Calling sequence ; CALL BALANC(NM,N,A,LOW,IGH, SCALE) 

Input arguments ; 

NM First dimension of array A as given in the calling program 

N Number of rows of matrix A 

A Matrix which is to be balanced stored in a real two-dimensional 

array. The contents of this array are destroyed upon return. 

Output arguments ; 

Upon normal return, A contains the balanced matrix. 

Two integers such that A(I,J) = 0 if I > J and 
J = 1 ,2, . . . ,L0W-1 or I = IGH+1 ,IGH+2, . . . ,N. 

A one-dimensional array of dimension at least N. SCALE contains 
information determining the permutations and scaling factors 
used; suppose that the principal submatrix in rows LOW through 
IGH has been balanced, that P(J) denotes the index inter- 
changed with J during the permutation step, and that the ele- 
ments of the diagonal matrix used are denoted by D(I,J); then, 

SCALE(J) □ P(J) (J = 1 ,2, . . . ,L0W-1 ) 

= D(J,J) (J = L0W,L0W+1 , . . . ,IGH) 

= P(J) (J = IGH+1, IGH+2 N) 

The order in which the interchanges are made is N to 
IGH + 1, then 1 to LOW-1. 

COMMON blocks ; None 

Error message s; None directly from BALANC 
Field length ; 327 octal words (215 decimal) 

Subroutines employed b y BAL ANC; None 
Subroutine em plo ying BALAN C : EIGEN 

Comments; None 


A 

LOW, IGH 
SCALE 
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Subroutine ELMHES 


Description ; Given a real square matrix A, the purpose of ELMHES is to 
reduce a submatrix situated in rows and columns LOW through IGH to upper 
Hessenberg form by stabilized elementary similarity transformations. The 
computational method follows that of Martin and Wilkinson in reference 18. 
ELMHES is a translation of the ALGOL procedure ELMHES found on pages 339 
to 358 of reference 2. 

Source of software ; LaRC Analysis and Computation Division subprogram library 
Calling sequence ; CALL ELMHES (NM,N, LOW, IGH, A, INT) 

Input arguments ; 

NM First dimension of the array A as given in the calling program 

N Number of rows of matrix A 

LOW, IGH Integers typically determined by the balancing subroutine BALANC 
for eigenvalue and eigenvector computation. If BALANC is not 
used, set LOW = 1 and IGH = N. 

A Matrix which is to be used in the reduction to upper Hessenberg 

form stored in a real two-dimensional array. The contents 
of A are destroyed upon return. 

Output arguments ; 

A Upon normal return, A contains the upper Hessenberg matrix. 

The multipliers which were used in the reduction are stored in 
the remaining triangle under the Hessenberg matrix. 

INT One-dimensional integer array of dimension at least IGH in the 

calling program. INT contains information on the rows and 
columns interchanged in the reduction. Only elements LOW to 
IGH are used. 

COMMON blocks ; None 

Error messages ; None directly from ELMHES 
Field length ; 206 octal words (13^ decimal) 

Sub routines employed by ELMHES ; None 
Sub routine employing ELMHES ; EIGEN 
Comments; None 
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Subroutine HQR 


Description ; The purpose of HQR is to find the eigenvalues of a real square 
upper Hessen berg matrix H by the QR algorithm. The computational method 
follows that of Martin, Peters, and Wilkinson found in reference 19- The 
subroutine HQR is a translation of the ALGOL procedure HQR found on pages 359 
to 371 of reference 2. 


Source of software : LaRC Analysis and Computation Division subprogram library 

Calling sequence : CALL HQR(NM,N,LOW,IGH,H,WR,WI,IERR) 


Input arguments : 


NM 


First dimension of array H as given in the calling program 


N 


Number of rows in matrix H 


LOW, IGH Integers typically determined by the balancing subroutine BALANC. 
If BALANC is not used, set LOW = 1 and IGH = N. 

H Matrix in upper Hessenberg form stored in a real two-dimensional 

array. Information about the transformations used in the reduc- 
tion to Hessenberg form by subroutine ELMHES, if performed, is 
stored in the remaining triangle under the Hessenberg matrix. 
Upon normal return, H is destroyed. 

Output arguments : 


WR, WI 


lERR 


One-dimensional arrays, dimensioned at least N in the calling 
program, containing, upon normal return, the real and imaginary 
parts of the eigenvalues, respectively. The eigenvalues are 
unordered except that the complex conjugate pairs of values 
appear consecutively with the eigenvalue having positive 
imaginary part first. 


Integer error code: 

lERR = 0 Normal return 

lERR = J The Jth eigenvalue has not been determined after 

30 iterations of the QR algorithm. If an error 
exit is made, the eigenvalues should be correct 
for indices IERR+1 ,IERR+2 , . . ..,N. 
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COMMON blocks : None 


Error messages ; None directly from HQR; lERR should be examined upon return. 
Field length ; 556 octal words (366 decimal) 

Subroutines employed by HQR ; None 
Subroutine employing HQR : EIGEN 


Comments: None 



Subroutine INVIT 


Description : The purpose of INVIT is to find those eigenvectors of a real 

square upper Hessenberg matrix corresponding to specified eigenvalues using 

inverse iteration. The subroutine INVIT is a translation of the ALGOL pro- 
cedure INVIT found on pages 41 8 to 439 of reference 2. 

Source of softwar e: LaRC Analysis and Computation Division subprogram library 

Calling sequence : CALL INVIT(NM,N,A,WR,WI,SELECT,MM,M,Z,IERR,RM1 ,RV1 ,RV2) 

Input arguments : 

NM First dimension of the array A as given in the calling 

program 

N Number of rows in A 

A Matrix in upper Hessenberg form stored in real two- 

dimensional array. Upon return, A is unaltered. 

WR, WI One-dimensional arrays dimensioned at least N by the 

calling program. WR and WI contain the real and 
imaginary parts, respectively, of the eigenvalues of 
the matrix. The eigenvalues must be stored in a manner 
identical to that of subroutine HQR. WI is unaltered 
upon return. WR may be altered since close eigenvalues 
are perturbed slightly in searching for independent 
eigenvectors . 

SELECT One-dimensional array of logical variables, dimensioned 

at least N by the calling program. SELECT specifies 
the eigenvectors to be found. The eigenvector corre- 
sponding to the Jth eigenvalue is specified by setting 
SELECT(J) to TRUE. Upon return, SELECT may be altered. 

If the elements corresponding to a pair of conjugate com- 
plex eigenvalues were each initially set TRUE, the program 
resets the second of the two elements to FALSE. 

MM MM should be set to an upper bound for the number of columns 

required to store the eigenvectors to be found. Note that 
two columns are required to store the eigenvector corre- 
sponding to a complex eigenvalue. 

RM1 , RV1 , RV2 Temporary storage arrays dimensioned at least N x N, N, 

and N, respectively, by the calling program 


Output arguments: 


M 


The actual number of columns used to store the eigenvectors 



z 


lERR 


Matrix dimensioned at least NM x MM by the calling program. 
Z contains the real and imaginary parts of the eigen- 
vectors. If the next selected eigenvalue is real, the 
next column of Z contains its eigenvector. If the eigen- 
value is complex, the next two columns of Z contain the 
real and imaginary parts of its eigenvector. The eigen- 
vectors are normalized so that the component of largest 
magnitude is unity. Any vector which fails the acceptance 
test is set to zero. 


Integer error code: 
lERR = 0 

lERR = -(2N + 1) 


lERR = -K 
lERR = -(N + K) 


Normal return 

More than MM columns of Z are 
necessary to store the eigenvectors 
corresponding to the specified 
eigenvalues. 

The iteration corresponding to the 
Kth value fails. 

Both of the above error situations 
occur. 


COMMON blocks: None 


Error messages : None directly from INVIT; lERR should be examined upon return. 

Field length : 1310 octal words (712 decimal) 

Subroutines employed by INVIT : None 

Subroutine employing INVIT : EIGEN 


Comments : None 
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Subroutine ELMBAK 


Description ; The purpose of ELMBAK is to form the eigenvectors of a real 

square matrix A by back transforming those of the corresponding upper 

Hessenberg matrix determined by subroutine ElUHES. ELMBAK is a translation 

of the ALGOL procedure ELMBAK found on pages 339 to 358 of reference 2. 

Source of software ; LaRC Analysis and Computation Division subprogram library 
Calling sequence ; CALL ELMBAK(NM,LOW,IGH,A,INT,M,Z) 

Input arguments ; 

NM First dimension of the array A as given in the calling program 

LOW, IGH Integers available from the balancing subroutine BALANC. If 

BALANC is not used, set LOW = 1 and IGH = The order of the 
matrix. 

A Two-dimensional array dimensioned at least NM x IGH in the 

calling program. A contains the multipliers which were used 
in the reduction by ELMHES in its lower triangle below the 
subdiagonal. 

INT One-dimensional array dimensioned at least IGH by the calling 

program. INT contains information on the rows and columns 
interchanged in the reduction by ELMHES. Only elements LOW 
through IGH are used. 

M The number of columns of Z to be back transformed 

Z Two-dimensional array dimensioned at least NM x M by the 

calling program. Z contains the real and imaginary parts 
of the eigenvectors to be back transformed in its first 
M columns. The contents of Z are destroyed upon return. 


Output argument ; 

Z The real and imaginary parts of the transformed eigenvectors are 

returned in the first M columns. 

COMMON blocks ; None 

Error messages ; None directly from ELMBAK 
Field length ; 133 octal words (91 decimal) 

Subroutines employed by ELMBAK : None 

Subroutine employing ELMBA K : EIGEN 

Comments; None 
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Subroutine BALBAK 


Description ; The purpose of BALBAK is to form the eigenvectors of a real 

square matrix by back transforming those of the corresponding balanced matrix 
determined by subroutine BALANC, BALBAK is a translation of the ALGOL pro- 
cedure BALBAK found on pages 315 to 326 of reference 2. 


Source of software ; LaRC Analysis and Computation Division subprogram library 
Calling sequence ; CALL BALBAK(NM,N, LOW, IGH, SCALE, M,Z) 


Input arguments ; 

NM First dimension of the matrix array as given in the calling 

program 

N Number of rows of the matrix 


LOW, IGH Integers determined by BALANC 

SCALE One-dimensional array dimensioned at least N by the calling 

program. SCALE contains information determining the permuta- 
tions and scaling factors used by BALANC. 


M 


The number of columns of Z to be back transformed 


Z Two-dimensional array dimensioned at least NM x M in the call- 

ing program. Z contains the real and imaginary parts of the 
eigenvectors to be back transformed in its first M columns. 


Output argument ; 

Z The real and imaginary parts of the transformed eigenvectors are 

returned in the first M columns. 

CO MMON blocks ; None 

Error messages ; None directly from BALBAK 
Field length ; 12M octal words (84 decimal) 

Subroutines employed by BALBAK ; None 
Subroutine employing BALBAK ; EIGEN 
Comments: None 
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Subroutine DETFAC 


Description ; The purpose of DETFAC is to factor a real square matrix A as 
PA □ LU 

where P is a permutation matrix representing row pivotal strategy, L is 
a unit lower triangular matrix, and U is an upper triangular matrix. 

Options are provided to compute the determinant of A with and without 
A input in factored form. 

Source of software ; LaRC Analysis and Computation Division subprogram library 
with modifications by Ernest S. Armstrong, LaRC 

Calling sequence ; CALL DETFAC (NMAX,N, A, IPIVOT.IDET, DETERM, ISCALE,WK,IERR) 

Input arguments ; 

NMAX First dimension of the array A as given in the calling 

program 

N Number of rows of matrix A 

A Two-dimensional array dimensioned at least NMAX x N in 

the calling program. A is the matrix to be factored. 

If the factored form of A is input, 

A = (L\U) 

should be used neglecting the unity elements of L and 
the pivotal strategy input through the array IPIVOT. 

For A in unfactored form, input data are destroyed. 

IPIVOT One-dimensional array dimensioned at least N by the call- 

ing program. Not required as input if the unfactored 
form of A is used. Otherwise, IPIVOT(I) = J indicates 
that row J of matrix A was used to pivot for the 
Ith column. 

IDET Determinant evaluation code ; 

0 Compute L and U matrices only. 

1 Given L and U matrices as input, compute param- 

eters defining the determinant. 

2 Compute L, U, and determinant parameters. 

WK One-dimensional array dimensioned at least N by the 

calling program and used as a work storage array. 
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Output arguments ; 


A Upon normal return, the L and U matrices are over 

stored in A as 

A = (L\U) 

neglecting the unity elements in L. 

IPIVOT Upon normal return, IPIVOT contains the pivotal strategy 

as previously explained. 

DETERM, ISCALE Determinant evaluation parameters; upon return, for 

IDET i 0, 

det(A) = DETERM x lO^OO^ISCALE 


lERR Singularity test parameter; 

0 Matrix A is singular. 

1 Matrix A is nonsingular. 


COMMON blocks ; None 

Error messages ; None directly from DETFAC; lERR should be examined upon return. 
Field length ; 332 octal words (218 decimal) 

Subrouti nes employed by DETFAC ; None 
Subrou tine employing DETFAC ; GELIM 
Comments ; None 
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Subroutine AXPXB 


Description ; The purpose of AXPXB is to solve the real matrix equation 
AX + XB = C 

where A, B, and C are constant matrices of order m x m, n x n, and 
m X n, respectively. The matrices A and B are transformed into real 
lower and upper Schur form (ref. 5), and the transformed system is solved by 
back substitution. The option is provided to input the Schur forms directly 
and bypass the Schur decomposition. 

Source of software : Coded directly from reference 5 

Calling sequence ; CALL AXPXB(A,U,M,NA,NU,B,V,N,NB,NV,C,NC,EPSA,EPSB,FAIL) 

Input arguments ; 

A Two-dimensional array dimensioned at least (m + 1) x (m + 1) 

in the calling program. The upper m x m part contains the 
matrix A. If the Schur form of A is input directly the lower 
triangle and superdiagonal of the upper m x m part of the 
array A contain a lower Schur form of A. 

U Two-dimensional array dimensioned at least m x m in the calling 

program. Not required as input if the Schur form of A is not 
input directly. Otherwise, U contains the orthogonal matrix 



that reduces 

A 

to Schur 

form 

A 

via 


A = 

U' 

AU 




M 

Number of rows 

of 

the matrix 

A 



NA 

First dimension 

of 

the array 

A; 

at 

least m + 1 

NU 

First dimension 

of 

the array 

u; 

at 

least m 


B Two-dimensional array dimensioned at least (n + 1) x (n + 1) 

the calling program. The upper n x n part contains the 
matrix B. If the Schur form of B is input directly, the 
upper triangle and subdiagonal of the upper n x n part of 
array B contain an upper real Schur form of B. 

V Two-dimensional array dimensioned at least n x n in the calling 

program. Not required as input if the Schur form of A is not 
input directly. Otherwise, V contains the orthogonal matrix 
that reduces B to Schur form B via 

B = V*BV 

N Number of rows of the matrix B 


in 

the 
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NB First dimension of the array B; at least n + 1 

NV First dimension of the array V; at least n 

C Two-dimensional array dimensioned at least m x n by the calling 

program; destroyed upon return 

NC First dimension of the array C; at least m 

EPSA Convergence criterion for the reduction of A to Schur form. 

EPSA should be set slightly smaller than 10~ ^ where is 

the number of significant digits in the elements of A. Set EPSA 
negative if a Schur form for A and transformation matrix U 
are directly input. 

EPSB Convergence criterion for the reduction of B to Schur form. EPSB 

should be set slightly smaller than 10~ ^ where is the 

number of significant digits in the elements of B. Set EPSB 
negative if a Schur form for B and transformation matrix V 
are directly input. 


Output arguments ; 


A 

Upon normal return, 

the 

array 

A 

contains 

the 

lower Schur 

form 

of 


the matrix A. 









U 

Upon normal return, 

the 

array 

U 

contains 

the 

orthogonal 

matrix 

U 


which transforms 

A 

to Schur : 

form. 





B 

Upon normal return. 

the 

array 

B 

contains 

the 

upper Schur 

form 

for 


the matrix B. 









V 

Upon normal return. 

the 

array 

V 

contains 

the 

orthogonal 

matrix 

V 


which transf orms 

B 

to upper ; 

Schur form. 





C 

Upon normal return. 

the 

solution 

X is stored 

in C. 




FAIL Integer variable containing an error signal. If FAIL is positive 

(negative) then the program was unable to reduce A(B) to real 
Schur form. If FAIL = 0, the reductions proceeded without 
mishap . 

COMMON blocks : None 

Error messages : None directly from AXPXB; FAIL should be tested upon return. 

Field length ; 66? octal words (439 decimal) 
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Subroutines employed by AXPlffl: HSHLDR, BCKMLT, SCHUR, SHRSLV 

Subroutine employing AXPXB : BARSTW 

Comments ; Subroutine AXPXB should be used when AX -i- XB = C needs to be 
solved for a number of C matrices. For the first C, set EPSA and EPSB 
based on kg^ and k^. Afterward, assuming PAIL = 0, set EPSA eind EPSB to 
negative values and compute X for the remaining C matrices. For the 
special case in which 

A = B» and C = C 

the subroutine ATXPXA should be employed. 
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Subroutine SHRSLV 

Description ; The purpose of SHRSLV is to solve the real matrix equation 
AX + XB = C 

where A is an m x m matrix in lower real Schur form and B is an 

n X n matrix in upper real Schur form. 

Source of software ; Coded directly from reference 5 
Calling sequence ; CALL SHRSLV(A,B,C,M,N,NA,NB,NC) 

Input arguments ; 

A Two-dimensional array dimensioned at least m x m in the calling 

program. A contains the lower Schur form of the matrix A. 

Not destroyed upon return. 

B Two-dimensional array dimensioned at least n x n in the calling 

program. B contains the matrix B in upper Schur form. Not 
destroyed upon return. 

C Two-dimensional array dimensioned at least m x n by the calling 

program. C contains the C matrix of the algebraic equation 
which is destroyed upon return. 

M Number of rows of the matrix A 

N Number of rows of the matrix B 


NA 

First 

dimension 

of 

the 

array 

A; 

at 

least 

m 

NB 

First 

dimension 

of 

the 

array 

B; 

at 

least 

n 

NC 

First 

dimension 

of 

the 

array 

c; 

at 

least 

ra 


Output argument ; 

C Upon normal return, the solution X is in C. 

COMMON block ; SLVBLK 

Error messages ; None directly from SHRSLV 
Field length ; MMM octal words (292 decimal) 

Subroutine employed by SHRSLV ; SYSSLV 
Subroutine employing SHRSLV ; AXPXB 
Comments ; None 
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Subroutine ATXPXA 


Description ; The purpose of ATXPXA is to solve the real matrix equation 
A'X + XA = C 

where A and C are constant matrices of dimension n x n with 


C = C‘ 


The matrix A is transformed into upper Schur form (ref. 5) and the trans- 
formed system is solved by back substitution. The option is provided to 
input the Schur form directly and bypass the Schur decomposition. 


Source of software ; Coded directly from reference 5 
Calling sequence ; CALL ATXPXAC A,U,C,N,NA,NU,NC,EPS,FAIL) 


Input arguments ; 


A 


U 


C 

N 

NA 

NU 

NC 


Two-dimensional array dimensioned at least (n + 1 ) x (n + 1 ) 
in the calling program. The upper n x n part contains the 
matrix A. If the Schur form of A is input directly the upper 
triangle and the first subdiagonal of the upper n x n part of 
the array A contain an upper real Schur form of A. 


Two-dimensional array dimensioned at least n x n in the calling 
program. Not required as input if the Schur form of A is not 
input directly. Otherwise, U contains the orthogonal matrix 
that reduces A to Schur form A via 

A = U'AU 

Two-dimensional array dimensioned at least n x n by the calling 
program; destroyed upon return 

Number of rows of the matrix A 


First dimension of the 
First dimension of the 
First dimension of the 


array 

A; 

at 

least 

n + 1 

array 

U; 

at 

least 

n 

array 

c; 

at 

least 

n 


EPS Convergence criterion for the reduction of A to Schur form. 

EPS should be set slightly smaller than 10” ^ where is 

the number of significant digits in the elements of A. Set EPS 
negative if a Schur form for A and transformation matrix U 
are directly input. 
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Out put arguments ; 


A Upon normal return, the array A contains the upper Schur form of 

the matrix A. Same as input if EPS < 0. 

U Upon normal return, the array U contains the orthogonal matrix U 

which transforms A to Schur form. Same as input if EPS < 0. 

C Upon normal return, the solution X is stored in C. 

FAIL Integer variable containing an error signal. If FAIL is nonzero, 

the program was unable to reduce A to real Schur form. 

If FAIL = 0, the reduction proceeded without mishap. 

COMMON blocks : None 

Error messages ; None directly from ATXPXA. FAIL should be tested upon return. 

Field length ; 562 octal words (370 decimal) 

Su broutine s empl oyed by ATXPXA : HSHLDR, BCKMLT, SCHUR, SYMSLV 

Subroutine employing ATXPXA ; BARSTW 

Comments ; Subroutine ATXPXA should be used when 
A'X + XA = C 

needs to be solved for a number of symmetric C matrices. For the first C, 
set EPS based on k^. Afterward, assuming FAIL = 0, set EPS negative and 
compute X for the remaining C matrices. 
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Subroutine SYMSLV 


Description ; The purpose of SYMSLV is to solve the real matrix equation 
A'X + XA = C 

where C = C* and A is n x n and in upper real Schur form. 

Source of software : Coded directly from reference 5 

Calling sequence ; CALL SYMSLV(A,C,N,NA,NC) 

Input arguments ; 

A Two-dimensional array dimensioned at least n x n in the calling 

program, A contains the upper Schur form for the matrix A. 
Not destroyed upon return. 

C Two-dimensional array dimensioned at least n x n in the calling 

program, C contains the C matrix of the algebraic equation 
which is destroyed upon return. 

N Number of rows of the matrix A 

NA First dimension of the array A; at least n 

NC First dimension of the array C; at least n 

Output argument ; 

C Upon normal return, the solution X is stored in C. 

COMMON block ; SLVBLK 

Error messages ; None directly from SYMSLV 
Field length ; 535 octal words (3^9 decimal) 

Subroutine employed by SYMSLV ; SYSSLV 
Subroutine employing SYMSLV; ATXPXA 


Comments : None 
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Subroutine HSHLDR 


Description ; The purpose of HSHLDR is to reduce a real n x n matrix A to 
upper Hessenberg form by Householder’s method of elementary Hermitian trans- 
formations described in reference 29. 

Source of software ; Coded directly from reference 5 

Calling sequence ; CALL HSHLDR(A,N,NA) 

Input arguments ; 

A Two-dimensional array dimensioned at least (n + 1 ) x (n + 1 ) in 

the calling program. The upper n x n part of the array A 
contains the matrix A which is destroyed upon return. 

N Number of rows of matrix A 

NA First dimension of the array A; at least n + 1 

Output argument ; 

A Upon normal return, the upper triangle of the array A to 

column n contains the upper triangle of the Hessenberg form 
of the matrix A. Column n + 1 contains the subdiagonal ele- 
ments of the Hessenberg form. The lower triangle and the 
(n + 1)th row of the array contain a history of the Householder 
transformations . 

COMMON blocks ; None 

Error messages ; None directly from HSHLDR 
Field length ; 265 octal words (181 decimal) 

Subroutines employed by HSHLDR ; None 
Subrout ines employing HSHLDR ; AXPXB , ATXPXA 
Comments ; None 
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Subroutine BCKMLT 


Description ; Given the output A from subroutine HSHLDR, the purpose of 
BCKMLT is to compute the orthogonal matrix that reduces A to upper 
Hessenberg form. 

Source of software ; Coded directly from reference 5 
Calling sequence ; CALL BCKMLT(A,U,N,NA,NU) 

Input arguments ; 

A Two-dimensional array containing the output from HSEILDR 

N Number of rows in matrix A in HSHLDR 

NA First dimension of the array A; at least N + 1 

NU First dimension of the array U; at least N 

Output argument ; 

U Two-dimensional array dimensioned at least n x n where n is 

the order of the A matrix in HSHLDR. If the matrix A is 
used for U in the calling sequence, the elements of the 
orthogonal matrix will overwrite the output of HSHLDR. 


COMMON blocks ; None 


Error messages ; None directly from BCKMLT 
Field length ; 201 octal words (129 decimal) 
Subroutines employed by BCKMLT ; None 
Subroutines employing BCKMLT : AXPXB, ATXPXA 

Comments : None 
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Subroutine SCHUR 


Desoription ; The purpose of SCHUR is to reduce an n x n upper Hessenberg 
matrix H to real Schur form. Computation is by the QR algorithm with 
implicit origin shifts. The program SCHUR is an adaptation of the ALGOL 
program HQR found in reference 19. The product of the transformations used 
in the reduction is accumulated. 

Sourc e of software ; Coded directly from reference 5 

Calling sequence ; CALL SCHUR (H,U,NN,NH,NU, EPS, FAIL) 

Input arguments ; 

H Two-dimensional array dimensioned at least n x n in the calling 

program. The array H contains the matrix H in upper 
Hessenberg form. The elements below the third subdiagonal are 
undisturbed upon return. 

U Two-dimensional array dimensioned at least n x n by the calling 

program. On input, U contains any square matrix desired. 

NN The number of rows in the matrices H and U 

NH First dimension of the array H; at least n 

NU First dimension of the array U; at least n 

EPS Number used in determining when an element of H is negligible. 

The element H(i,j) is negligible if |H(i,j)| S EPS x (|H|| 
where 1|H|| denotes the norm of H. 

Output arguments : 

H Upon normal return, H contains an upper Schur form of H. 

U Upon normal return, U contains the product of the input U 

right multiplied by the accumulated orthogonal transformations 
used to reduce H to Schur form. If the identity matrix is 
input for U, then the Schur form is 

U'HU 

FAIL Integer variable containing an error signal. If FAIL is positive, 

then the program failed to make the (FAIL-1) or (FAIL-2) 
subdiagonal element negligible after 30 iterations of the 
QR algorithm. 
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COMMON blocks : None 


Error messages : None directly from SCHUR 

Field length ; 520 octal words (336 decimal) 
Subroutines employed by SCHUR : None 

Subroutines empl oyi ng SCHUR ; AXPXB, ATXPXA 
Comments: None 
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Subroutine SYSSLV 


Description ; The purpose of SYSSLV is to solve the linear system 
Ax = b 

where A is an nxn(nS5) matrix and' b is n-dimensional vector. 
Solution is by Grout reduction. The matrix A, the vector b, and order n 
are contained in the arrays A, B,; and the variable N of the COMMON block 
SLVBLK. The solution is returned in the array B. 

Source of software ; Coded directly from reference 5 

Calling sequence ; CALL SYSSLV 

Input arguments ; None 

Output arguments ; None 

COMMON block ; SLVBLK 

Error messages ; None directly from SYSSLV 
Field length ; 227 octal words (151 decimal) 

Subroutines employed by SYSSLV ; None 
Subroutines employing SYSSLV ; SHRSLV, SYMSLV 
Comments ; None 
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Subroutine GAUSEL 


Description : The purpose of GAUSEL is to solve a set of linear equations, 


AX = B 


by the method of Gaussian elimination. The constant matrices A and B are 
of dimension n x n and n x r, respectively. No information is returned 
on the pivotal strategy or the value of the determinant of A. 

Source of software : LaRC Analysis and Computation Division subprogram library 

Calling sequence : CALL GAUSEL(MAX,N,A,NR,B,IERR) 

Input arguments : 

MAX First dimension of the array B given in the calling program 

N Number of rows of the matrix A 

A Two-dimensional array dimensioned at least n x n in the calling 

program; destroyed upon return 

NR Number of columns in the matrix B 

B Two-dimensional array dimensioned at least MAX x NR in the 

calling program; destroyed upon return 

Output arguments : 

B Upon normal return, B contains the solution X. 

lERR Integer error code: 

0 Normal return 

2 Input matrix A is singular. 

COMMON blocks ; None 

Error messages ; None directly from GAUSEL. lERR should be checked upon 
return. 

Field length : 317 octal words (207 decimal) 

Subroutines employed by GAUSEL : None 

Subroutine employing GAUSEL : EXPADE 

Comments: None 
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EXAMPLE COMPUTATIONS 


This section contains selected examples illustrating the use of the ORACLS 
program subroutines. Sample executive programs and output data are presented 
which develop state variable feedback control laws using the optimal transient 
regulator, optimal sampled-data regulator, and model-following design 
approaches. Additionally, the construction of an asymptotic Kalman-Bucy 
estimator is illustrated. 

Data employed for the plant equations are from a linearized mathematical 
model of the lateral dynamics of an F-8 aircraft, presented in reference 30. 

The problem and machine-dependent accuracy and convergence parameters required 
for COMMON blocks TOL and CONV of subroutine RDTITL used in the example compu- 
tations are: 

EPSAM = EPSBM = 1 .E-10 

lACM = 12 

SUMCV = 1.E-8 

MAXSUM = 50 

RICTCV = 1.E-8 

SERCV = 1.E-8 


The construction of RDTITL for the example computations is shown in figure 1 . 
All computations were performed using the Control Data Cyber digital computer 
system under Network Operating System (NOS) 1.2. 


SUBROUTINE RDTITL 
C 

COMMON/LINES/NLP, LIN, TITLE(8),TIL(2) 

COMMON/FORM/NEPR.FMTl (6) ,FMT2(6) 

COMMON/TOL/EPSAM.EPSBM, I ACM 
COMMON/CONV/SUMCV, MAXSUM, RICTCV.SERCV 
C NPL = NO. LINES/PAGE VARIES WITH THE INSTALLATION 
DATA LIN, NLP/1 ,44/ 

DATA NEPR,FMT1/7,10H(1P7E16.7)/ 

DATA TIL/10H 0RACL,10HS PROGRAM/ 

DATA FMT2/10H(3X,1P7E16,10H.7) / ' 

DATA EPSAM/l.E-10/ 

DATA EPSBM/1.E-10/ 

DATA IACM/12/ 

DATA SUMCV/1.E-8/ 

DATA RICTCV/l.E-8/ 

DATA SERCV/l.E-8/ 

DATA MAXSUM/50/ 

READ(5,100) TITLE 
IF(E0F(5))90,91 

90 CONTINUE 
STOP 1 

91 CONTINUE 
100 F0RMAT(8A10) 

CALL LNCNT(IOO) 

RETURN 

END 

Figure 1.- Subroutine RDTITL for the example computations. 
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Example 1 - Optimal Transient Regulator 


This problem illustrates the construction of storage arrays for ORACLS 
and demonstrates the solution of a simple transient optimal linear regulator 
problem. sGiven the system 


x(t) a A x(t) + B u(t) 


with 


and 


A = 


-2.6 

0.25 

-38. 

0.0 

-0.075 

-0.27 

4.4 

0.0 

0.078 

-0.99 

-0.23 

0.052 

1.0 

0.078 

0.0 

0.0 



17. 

7.0 


0.82 

-3.2 

B = 

0.0 

0.046 


0.0 

0.0 

control law 

u(t) = -1 

1 

J = - x»(20) 

x(20) + 


j; 


20 


The solution can be found by directly employing the subroutine CNTNREG. Output 
data are to be printed at 2.0-second intervals between 0 and 20. 

Referring to the description of CNTNREG, the dimensions of the subroutine 
argument arrays must be specified. The matrices A, B, Q, and R are 4x4, 
4 X 2, 4 X 4, and 2x2, respectively. As packed one-dimensional arrays, they 
must be dimensioned at least A(l6), B(8), Q(16), and R(4). However, for this 
example, input data to CNTNREG are defined within the source program and, for 
convenience of construction, two-dimensional arrays A(4,4), B(4,2), Q(4,4), 
and R(2,2) are used. Since H s I4, no data for H are required if the 
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logical variable IDENT is set to TRUE. The array P initially contains the 

1 

weighting matrix - I 4 on the final states and hence, is dimensioned 


The eirray F will contain the feedback gain matrix and must be dimensioned at. 
least F(8). Similarly, Z, W, LAMBDA, and S are dimensioned 64, 6^, 16, 
and 16 , respectively. Dimensions for NA, NB, NQ, NR, NF, and NP are set 
to 2. The vector T contains the final problem time and print interval and 
is also dimensioned 2. The vector lOP controls the printing from CNTNREG and 
must be dimensioned at least 3. The dimension of DUMMY is computed from the 
formula (since I0P(3) = 0 for a transient solution to the Riccati equation). 


9n^ + 17n + nr 


with n = 4 and r = 2. 

The following executive program shows how the specific data for CNTNREG 
may be constructed. Note that subroutines NULL and UNITY may be employed to 
generate zero and unity elements in the coefficient and weighting matrices. 
For the NOS 1.2 system at LaRC tape 5 is designated input and tape 6 output, 
as indicated by the source program card. Output from ORACLS follows the 
executive program. 





PfcnOSAM Rf^GlAT OPT»l 


FTN A.6+439 17 / 07 / 27 . 08.37.23 


1 


5 


10 


15 


30 


?5 


30 


35 


PP0GP4M REGlAK INPUT# OUTPUT. TAPS 5-INPUT# TAPE6«nUTPUn 

c 

OIMtN<;iaN A(4#4). 8(4#2)#0(4#4)#R(2#2)>Z{64)>W(64)#LAMB0Aa6)#S(l6) 
l#P(2#4>#P(16)#DUMMY(220).NA(2),N8(2)#NQ(2)#NR(2)#TOP(3)#T(2)#NF(2) 
2 # N P ( 2 ) 

LC’GICAL IDSMT 
RfAL LAMBDA 

C 

C 

r IMTULIZE NA.NP# •••#NP 

NA ( 1 4 
N'A(2>- 4 
NP ( 1 )« 4 
NP(2>- 2 
N0( 1 )« 4 
N0(2>* 4 
NP<1)« 2 
NM2>« 2 

MP<1). 4* 

NP«2)» 4 

C 

C SFT FINAL TIME AND PRINT INTERVAL 

T(l)-20. 

T(2) - 2. 

C 

C DEFINE PRINT AND TRANSIENT SOLUTION OPTIONS 

lOP(l) • 1 
irP(2> - 1 
IDP(3) • 0 

C 

C DEFINE COEFFICIENT AND WEIGHTING MATRICES 
CALL NULL (A,NAJ 
CAU NULUP.N3) 

CALL UNITY(P,NP) 

CALL SC ALE( P.NPf P#NP#0.5) 

CALL UNITY(O.NO) 

CALL UNITY(R.NR) 

CAIL SCALE<R#NR,R,NR#100. ) 
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C 


A(1,U« -2,6 
Afl,2) * .25 
A(l,3)»-38. 



PROGRAM REGLAT 


74/74 OPT-1 


PTN 4.6+439 


77/10/26. 18.59.15 



A(2,l) 

■ 

-.075 


A(2,2) 

m 

-.27 

4!? 

A12,3) 

M 

4.4 


A(3,l) 

m 

• 07b 


A( 3,2) 

s 

-.99 


A(3,3) 

3 

-.23 


A(3,4) 

3 

.052 


A < 4 , i ) 

3 

i.U 


A(4,2) 

3 

.076 


b(l,l) 

3 

17. 


B(l,2) 

3 

7.0 


B(2,l) 

3 

.82 

55 

B(2,2) 

3 

-3.2 


6(3^2) 

3 

• 046 


IDENT 

■ • 

TRUE. 


C 

C INPUT HOLLERITH DATA FOR TITLE OF CUTOUT 

60 CALL RDTITL 

C 

C NOW USE CNTNREG TO SOLVE THE TRANSIENT RFGULATOR PROBLEM 

CALL CNTNREG (AfNA#B#N6.H,NH,C, NO, P,NR, Z#W# LA MR D A> S# F # NF , P# NP, T# I OP 
1# IDENT,OUMMY) 

65 C 

C 

STOP 

END 


Jr 
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EXAMPLE X - OPTIMAL TRANSIENT REGULATOR 


QRACLS 


PROGRAM TO SOLVE THE T1 ME -INVAR lANT FINITE-DURATION CONTINUOUS OPTIMAL 
REGULATOR PROBLEM WITH NOISE-FREt MEASUREMENTS 


A MATRI 
-2*60000t0E+u0 
-7.500C000E-02 
7.8000000E-02 
1.0000000E400 


A ROWS 

2.SOOJLOC6-U1 

-2#7UOOCOOE-Oi 

-<i,900o00CE-Ol 

7.3000000E-02 


4 COLUMNS 
-3.3000000t+01 
4«40000N^0t 
-2.3000000E-CX 


0 . 

0 • 

5.2COOOOOE-02 

C. 


B matrix • 4 ROWS 2 COLUMNS 

1.7OCCOC0E+OX 7*0000t0u£+0v 
8.2000000E-01 -3,2C0JU00E+00 

0* 4.6CCOOOOE-U2 

0 « 0 . 


Q MATRIX 

4 ROWS 

4 COLUMNS 

l.OOOOOOOf+OJ 

u « 

0* 

0. 

1. JbOJOOOt + Ou 

u • 

0. 

0. 

i.OOCOOOuE+OC 

0« 

0. 

0. 


0 * 

0 * 

0 • 

l.OOOOOOOE+00 


H IS AN IDENTITY MATRIX 


R MATRIX 2 ROWS 2 COLUMNS 

1.0C00000E+U2 0* 

0 . X.0000u0CE4U2 


WEIGHTING ON TERMINAL VALUE OF STATE VECTOR 


P MATRIX 

4 ROWS 

4 COLUMNS 


5.0000000E-01 

0« 


c« 

c. 

5*OCOOOOOE-Ol 



0. 

0. 

5,0000000E-01 

0« 

0« 

0* 


5.OOG0OCOE-O1 


PROGRAM 


Z 


matrix 


8 ROWS 


6 COLUMNS 



EXAMPLfc 1 - OPTICAL TRANSIENT REGULATOR 


ORACLS PRQGRAH 


-2t 6000000 E+00 2 . &OOOOOOE-01 

0* 

-7.5000000E-02 -2 . 7uOOOOOE-Ol 

0. 

7.8000000E-02 -9. 900J0O0E-01 

0* 

l.OOOOOOOt+00 7*30O0C0CE-O2 

0 • 

^l.OOOOOCOEfOO c* 
*1«OOOOOOOE+Oo 

0. -l.OOCOCOOE+OO 

-7.6000000E-02 

0. 0« 

0* 

0* 0. 

0. 


8«A600000€-02 

*I«0912400£-01 

1«4720000E-03' 

0. 

7.5000000E-02 

2.7000000E-01 

^^•AOOOOOOE^OO 

0» 


*a»2200000E-03 

1.^7200008-03 

-2.1160000E-05 

0. 

-7.8000000E-02 

9.9000000E-01 

2.3000000E-01 

-5.2000000E-02 


EIGENVALUES OF I 


-6.4605654E-01 

6.4605654E-01 

7.8466995E-01 

7.8466995E-01 

-7.8466995E-U1 

-7.8466995E-01 

-2.8451457E+0O 

2.8451457E+00 


0. 

0. 

2.6501143E+00 
-2.65C1143E+00 
2.650U43E + 00 
-2.6501143E+00 

0. 

0. 


-3.80C0000E+01 

4.4OC000OE+00 

-2.30000aOE-01 

0. 

0« 

0 . 

-l.UOOOOOOt+OO 

0. 


!>.0257395E-02 

-1.634A167L-02 

1.7516130E-06 

1.^5U446E-C2 


0. 

0. 

5.2C00000E-02 

0. 

0. 

0. 

0. 

-I.OCOOOOOE+00 


2.6582832E-02 

9.5183892E-03 

-8.0834088E-03 

-1.4186503E-02 


-3.3800000E+00 

3.4600000E-02 

-3.2200000E-03 

0 . 

2.60000008^00 

-2.5000000E-C1 

3.80000008401 

0 . 


-8.5528327E-02 

3.5983643E-02 

-4.9824656E-03 

-4.8131275E-02 


-1.6423055E-01 

1.2843069E-02 

1.7743854E-02 

4.S46949SE-02 


6.2134586E-01 

3.0404804E-02 

-2.6968169E-03 

-2.1922162E-01 


CORRESPONDING EIGENVECTORS 


-2.4733021E-01 -1. 7233188E-01 

-6.2618216E-02 

4.3326566E-03 -S.3791761E-02 

-6.9930617E-03 

7.5286043E-03 7. 7409397E-03 

3.3935262E-04 

3.8230750E-01 -2 . 7082403E-0i 

-2.22005C7E-02 



EXAMPLl 1 - OPTICAL TRANSIENT REGULATOR 


QRACLS PROGRAH 


7. 255858 AE-02 9. 5 50 7170E -02 

1.0AA1627E-01 

5.3890433E-ai 7 • 9262267fc-01 

3. A98253AE--J1 

-A.32O7O09E-O1 3.2205187E-01 

9,28532A5c-0i 

5.5697883E-0i 3.932 7A1AE -01 

-9.167608OE-03 


REORDERED EIGENVECTORS 


-X. 72331 0 8E-Oi 5.02 5 739 51.-02 

6.213A566E-01 

-3.379l76iE-u2 -1 . 63A^167E-02 

3.0A04804E-02 

7.7909397E-f3 1 • 7 5 16 130E-06 

-2.6988189E-03 

-2. 7082^03 E-Oi leA5UA^tt-02 

-2.1922162t-01 

9.5507170E-02 -2 . 2 37 76 AOE-02 

0. A269310E-02 

7.926 22 6 7E-01 3 . 3965 A 17 1-01 

2.37A9A20E-01 

3. 2205187E-01 J . 7A3 52 13E-01 
-7.0239213E-01 

3.9 327<»iAE-Qi -1 . <» 1206 36E -02 

-6.98885A0E-O2 


LAHBOA MATRIX OF EIGENVALUES OF 


0 . 


6.A605654E-01 

0 . 

0 . 

0 . 

WIZW MATRIX 


7.8A66995E-01 
-2.650ii<»3E + UU 

0 . 

6 ROWS 


-2.237764CE-02 

3.3965417E-01 

1.79352i3b-01 

-1.^120636E-C2 


2.6582032E-C2 

9.5183892E-C3 

-8.083A088E-03 

-i.<*ie6503E-02 

2.AA83836E-02 

-7.86336J1E-03 

9.2129547E-01 

A.7159165E-03 

Z v«ITH POSITIVE 

0 . 

2.65011A3E+00 

7.8A66995E-01 

0 . 

8 COLUMNS 


2.A483836E-02 

-7.8633601E-93 

9.21295^7E-Ol 

A.7159165E-03 


-6.2618218E-02 

-6.9930617E-03 

3.3935252E-04 

-2.2200507E-0? 

1.0A41627E-01 

3.4982534E-01 

9.2653245E-01 

-9.1676080E-03 

REAL PARTS 

0 . 

0 . 

0 . 

2.8451457E+00 


-1.5780566E-02 

3.2328081E-01 

-’•95O6394E-0X 

-3.8062689E-02 


-2.473302XE-01 

4.3326566E-03 

7.5286043E-03 

3.023O75OE-OX 

7.2558584B-02 

5.3890433E-OX 

-4.3207089E-01 

5.5697883E-0X 


-4.4315091E-02 

-7.4606776E^03 

8.7411594E-01 

-1.26B1699E-02 


-X.6423055E-01 

Xo2843065E-02 

X.7745854E-02 

4.5465498E-02 

-4.43X5091E-02 

-7.4606776E-03 

8#7411594E-01 

-1.268X699E-02 


8.42693X0E-O2 

2.3749420E-01 

-7.02392X3E-0X 

-8.98B8548E-02 


-8.5528327E-02 

3.5983643E-02 

-4.9824656E-03 

-4.8X31275E-02 

-X.5780566E-02 

3.232808XE*0X 

-2#9$06394E-0X 

-3.6062689E-02 
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6.4605654E-01 A.3S3909AE-13 -1.6913328E-13 

1.7A80680E-13 

9*9622886E-1<» 7. 8<»6699J>E-01 2.65011A3E + 00 

1.1224239E-13 

l»5597072£-i3 -2.650UA3E + 00 7.8A66995E-01 

-1.936966XE-12 

-6.5AA6563E-1A -9, 1129A^9E-13 -7. A776612E-1 3 

-1.2551633E-12 

6.0327566E-14 -1.2656352E-12 -3 . 2722 815E-13 

-1.A901340E-12 

1.3255405E-13 -3. 0IA7102E-12 -1 . A7A0753E-1 3 

-3.A2A7A65E-12 

-2.2830A83E-L3 1 . 322 7C 21E-12 3 . 0763 A66L-1 2 

9.5261A12E-13 

1.69A08 85E-13 - 1 . 769t>8 65 E-1 2 “5«J5335o3t-13 

•2.6A61A57E+JO 

Wll MATRIX A ROWS A COLUMNS 

-1.72331 68 E-01 5 . 025 7396E-02 2 . 656 2 8 32E-C2 

-3.3791761E-02 -1 . 63 A A167E-02 9 . 5 18 3892E -03 

7.7A09397E-U3 1 ♦ 75 I6l 30E-06 -8.083A038E-C3 

-2.7082AU3E-01 l.A5ilAA6E-o2 -1.A1665J3E-02 

W21 MATRIX A ROWS A COLUMNS 

9.5507170E-02 -2 . 2 3776 ACE -02 2 • A A 6 3 8 36E -02 

7.9262267E-01 3. 3965A17E-01 -7.8633601L-03 

3.2205187E-01 1.7A3:?213E-01 9.21295A7F-01 

3.9327A1AE-01 -1 . A12 J636E-02 A. 7i59l65b-03 

W12 MATRIX A ROwS A COLUMNS 

-2.A733021E-01 -1 .6A23055E-01 -8 . 5528327E-C2 

A.332 6566E-0 3 1 . 2 8 A 3065 E -U 2 3 . 5 98 3 6 A3E-C 2 

7.52860A3E-03 1. 77A565AE-02 -A. 982 A6p6t-C3 

3.62307506-01 A . 5 A65 A96E-02 -A . « 1 3 12 75 E-02 

W22 MATRIX A ROWS A COLUMNS 

7%255858AE-02 -A.A315L91t-02 -i . 57 80 566E-02 

5.3890A33E-01 -7. A606776E-03 3. 2 32 808 lE-01 

-A.3207089t-01 6 . 7A1159AE-01 -2.950639AL-C1 

5.5697883E-01 -1 . 2 681699E-02 -3 . 8 062 689E-02 

S MATRIX A ROWS A COLUMNS 


Jr 


QRACLS PROGRAM 


-2.0390208E-13 

-A.8159991E-13 

2.6A22617E-IZ 

2.8A51A57E+00 

1. A561670E-T2 

2. A210317E-12 
-1.7386355E-12 

2.5116A53E-1? 


-7.32A7951E-1A 

-8.27A6769E-13 

-1.7036632E-13 

-1.3339999E-X3 

-6.A60565AE-01 

-3.AA1A850E-X3 

1.0132695E-12 

-2.5555875E-13 


-3.A396XA0E-X3 

4.7260X36E-X3 

X.X402066E-X2 

4.7238794E-X4 

8.X668X50&7X3 

-7.8A66995E-0X 

2.650XXA3E4>00 

7.2909A086-X3 


8.6X90507E-X4 

X.5320X96E-X2 

-X.4882895E-X2 

6.XX8757XE-X3 

-4.39387I7E-X3 

-2«6$0XX43E^00 

-7#8466995E-0X 

-X*8733266E-12 


-6.2618216E-02 

-6.9930617E-03 

3.3935252E-0A 

-2.22C0507E-02 


1.0AA1627E-01 

3.A98253AE-01 

9.28532A5E-01 

-9.16760eOB-03 


6.213A566E-01 

3.0A0480AE-0? 

-2.6988i69E-03 

-2.1922162E-01 


8.A269310E-02 

2.37A9A20E-01 

-7.0239213E-01 

-8.98885488-92 



EXAHPLE I 


OPTINAL TRANSIENT REGUli^TOR 


ORACLS PROGRAH 


6t«>197403E-Ol 

9.301A743E-0X 

I«0292420£*00 

A.1882228E-01 


e«80<»0174E-02 

l«i451662E-0X 

-6.8625132E-01 

•2#9967598E-0X 


-/•06X7792E-03 

-^•2235758E-01 

8#0290849E-0X 

-4%0340X42fc-0X 


NATRIX (R 1NVERSE)X(6 TRANSPOSE) 


1.7000000E-OX 

7*0000000E--02 


EXP(-LANBDA X 


2.7468973E-OX 

0* 

0. 

0* 


6.200000CE-03 

-3.2000000E-02 


•20000000E401) 


0* 

1*X545067E-0X 

-X.7323710E-OX 

C. 


0* 

4«600Q0JCE^04 


0* 

X*73237XbE-01 

i,X545067fc-0I 

0. 


TIME ■ «20COOOO(rE-t>C2 


P MATRIX 
5*OOOOOOOE-OX 

0. 


0« 

0. 


4 ROWS 

0 • 

S.OOOOOOOE-OX 


0 . 

c. 


4 COLUMNS 

0. 

0 • 

5.000000CE-0X 

0. 


TIME ■ .X800000LE+02 


P MATRIX 
2.24806X1E-0X 
3.6054050E-0X 
-X.4217206E+00 
3.202386XE-0X 


4 ROWS 

3.60S40&CE-OX 

I.0744607E+Oi 

-6.86&5022E+00 

1#54C9X3SE+00 


4 COLUMNS 
-X.42X7206E+0C 
•a*8655022E+0C 
4.5683890t+01 
-2.6235U8fc + C0 


-1.0579571E-01 

-X.6582063E+00 

•1.3835163E-01 

1.2407188E+00 


0 . 

0. 


0 . 

0. 

0* 

3.3786086E-03 


0* 

0« 

0. 

5,0000000E-01 


3.202388XE-01 

X.5409X35E+00 

-2.6235XieE+00 

X#564XX09E+00 



EXAMPLE 1 - OPTIMAL TRANSIENT REGULATOR 


F MATRIX 2 ROWS COLUMNS 

4.1173471E-02 1.^939767E«01 -3 . 1438962E-0i 

3.5451406E-03 -3* 2266773L-01 2.0519022E-0X 


TIME • •16000000E4’02 


P MATRIX 4 ROWS 4 COLUMNS 

2.3819322E-01 4.0985673E-01 -1 • 5275096E+00 

4.0985673E-01 1»0991594E+01 -9. 0079421E+00 

-1.5275098E+00 -9. 007942 lE+00 4. 925632 8E+01 

3.6636522E-01 !• 713665 lE'fOO -2.9811626E'('C0 

F MATRIX 2 ROWS ^ COLUMNS 

4*3853673E-02 l*5980672fc-01 -3. 335418CE-01 

2.8554559E-03 -3#27X3471£-0l 2*0398637E-01 


TIME ■ .1400000CE+C2 


P MATRIX 4 ROWS 4 COLUMNS 

2.3924330E-01 4. 1302539E-01 •1.5297598E+00 

4#1302539E-0i !• 1028965E+01 -9.007650eE+00 

-1.5297598E+00 -9, 0076508E+00 4.9271952E+01 

3.6980853E-01 1. 7278^48E+00 -2 •98770596+00 

F MATRIX 2 ROWS 4 COLUMNS 

4#4058169E-02 1.6C65183E-01 -3* 3392190E-01 

2t8265289E-03 -3.2815862E-01 2.0362674E-01 


TIME • .12000000E+02 


P MATRIX 4 ROWS 4 COLUMNS 

2.3930520E-01 4 • 1337855E-01 -1 • 5297546E+00 




ORACLS PROGRAM 


6.7076089E-02 

-2*8099331E-02 


3.6636522E-01 

1.7136851E+00 

-2.9611826E+00 

1.7228142E+00 


7.6334306E-02 

-3.0563703E-02 


3*6980853E-0l 

1.7278448E+00 

2,9877059E+00 

1.7346173E+00 


7.7035778E-02 

•3*0778783E-02 


3.7003289E-01 



V71 
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EXAMPLE i - OPTIMAL TRANSIENT REGULATOR 


ORACLS 


A.l33765!>E-wl 1. 1O32206E + 01 

-1.52975A6E+00 -9.0C96623E400 

3.7003289E-01 i ,729320^6 +00 


-9#0096823E+00 

<r.92759UE+Cl 

-2.98e02i8E+00 


1.729320AE+00 

-2.9860218E+00 

1.735A617E+00 


F MATRIX 2 ROWS 

A.A071588E-02 1.0O7384AE-01 

2. 819563 lE-03 -3. 2623855E-01 


A COLUMNS 

-3.339376bt-Ol 7.7C86018E-02 

2, 036939 3E-01 -3.0810A39E-02 


TIME - .lOOOOOOCE+02 


P MATRIX 
2.39309X7E-01 
A*13A1056E-01 
-1.5297727E+00 
3.700A881E-01 


A ROWS 
A.13A1C56E-C1 
l,lC32A69E+0i 
-9.0098261E+00 
l,729AA93E+00 


A COLUMNS 
-1.5297727E+0C 
-9.0098261E+00 
A.9276275E+01 
-2.908O973E+OO 


3.700A881E-01 

1.729AA93E+00 

-2,9860973E+00 

1.7355256E+00 


F MATRIX 2 ROWS 

A.A072525E-02 i. 607A60AE-01 
2.8188087E-03 -3. 282AA78E-01 


A COLUMNS 

-3.339A193E-01 7. 7089781E-02 

2.03897A3E-C1 -3.0813A85E-02 


TIME ■ .80000000E+01 


P MATRIX 
2.39309A8E-U1 
A.13A1291E-01 
-1.52977A9E+00 
3.7005U05E-01 


A ROWS 
-*.13Ai29ie-0i 
1.1032A89E+01 
-9,OC93A19b+00 
1.729A59GE+U0 


A COLUMNS 
-i.52977A9t+00 
-9.0098A19E+00 
A.9276293E+01 
-2.9881061E+00 


3.7005005E-01 

1.729A590E+00 

“2.9881061E+00 

1.7355306E+00 


F MATRIX 
A.A072598E-02 
2.8X875A0E-03 


2 ROWS 

X.6C7A661E-01 

-3.282A52ttE-01 


A COLUMNS 
-3.339A2AAE-01 
2.0389779E-0i 


7.7090072E-02 

-3.081371AE-02 


PROGRAM 


TIME • 


600000 obE+OX 



example 1 - OPTIMAL TRANSIENT REGULATOR 


p 

MATRIX 

A ROWS 

^ COLUMNS 

2#3930951E-01 

4.1341309F-01 

-1.5297751E+00 

A,13<*1309E-01 

1.1032<»91E4‘01 

-9,0098<i33E + 00 

1.5297751E+00 

-9.0098<»33E+00 

<».927629AE + 01 

3*70050I5E-01 

1.729A597E+00 

-2,9081O68E+OO 

E 

MATRIX 

2 ROWS 

A COLUMNS 

«.407260<»E-02 

l,507^665E-0l 

-3.339A2^8E-01 

2.6187502E-03 

-3.282A531E-01 

2.0389783E-01 


ORACLS 


3.7005015E-01 

1,729^597E*00 

-2.9081O68E+OO 

1#73553I0E*00 


7.709009AE-02 

-3.0813729E-02 


.^OCOOOOOE+01 


TIME ■ 


P MATRIX 
2.3930951E-01 
4.13AX310e-Cl 
-1,5297751E^OO 
3.70050155-Cl 


F MATRIX 
^•^07260<»E-02 
2.8187499E-03 


^ ROWS 

^..13^1310F-01 
1.1032^i91F>01 
-9.009P<»3^E + 00 
1.729A597t+00 

2 ROWS 
1,607^665E-01 
-3,282<.531E-01 


COLUMNS 

-1.5297751E^00 
-9.0098A3«»E + 00 
<».927629AE + 01 
-2.9881068E+00 

^ COLUMNS 
-3.339^2<t8E-01 
2.0389783E-01 


3.7005015E-01 

l,729<^597E+00 

-2.9881068E+00 

l,7355310E+00 


7.7090096E-02 

-3.0813730E-02 


PROGRAM 


STEADY-STATE SOLUTION HAS BEEN REACHED IN CNTNREG 



Example 2 - Optimal Sampled-Data Regulator 

This problem demonstrates a situation in which several ORACLS subroutines 
are required to obtain a solution and input data are entered through use of the 
READ subroutine. 

Given the linear time- invariant system 


x(t) = A x(t) + B u(t) 


with x(0) □ xq given and A and B as defined in example 1, find the con- 
trol law 


u(t) = -F x(t) 


which minimizes 


J 


lim { r fx’Ct) Q x(t) + u*(t) R u(t)] 



where Q and R are 4x4 and 2x2 identity raalxrices, respectively, and 
u(t) is restricted to be piecewise constant over uniform time intervals of 
measure A = 0.5 second. Under these restrictions, the dynamical equations 
and performance index J become 


x[(i+1)Aj = A x(iA) + B u(iA) 


and 


J 


f ^ 

lim/ y [j 


x*(iA) Q x(iA) + x'(iA) W u(iA) + u*(iA) 



where 

A = e^ 


B = \ e*’’- B dT 

•^0 


i-52 



and 



Q dT 


Q H(t,0) dT 


H'(t,0) Q H(t,0 


B dT 


)] dT 


Performing the control variable transformation, 


u(iA) = -F x(iA) + v(iA) 


with 



2 


leads to 


x[(i+1)A] = A x(iA) + B v(iA) 


r N 


J = lim' 


a [x'(iA) Q x(iA) + v‘(iA) 


A = A - BF 


/V 

B = B 


where 


»> 




R = R 


Ignoring the initial value x(0) in J gives the problem of minimizing 

( N N 

[(i+1)A] Q x[(i+1)A3 + v*(iA) R v(iA)j 
i=0 , 


subject to 


x[(i+1)^ = A x(iA) + B v(iA) 
whose solution, if it exists, is given by 
v(iA) = -F x(iA) 

The gain matrix for the original problem is then 
F = F + F 


The complete problem can be solved by appropriately combining the subrou- 
tines^ EXPINT (for A and |) , SAMPL (for Q, W, and R) , PREFIL (for A, Q, 
and F) , and ASYMREG (for F) . An executive program for constructing the solu- 
tion follows. Data are input by using subroutine READ rather than being defined 
internally to the source program, as in example 1 . The dimension of DUMMY 
is chosen to be the maximum of the requirements of each operation employing 
this storage vector. At the stage where matrices F and F are printed, the 
statement ’’CALL LNCNT (100)” is used to shift the printing to the top of the 
next page. Output data are presented after the executive program. 
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PROGRAM SAMOAT 


74/74 OPT-1 


FTN 4.6+43<> 


77/07/27. 08.36.31 


X 


5 


10 


15 


20 


25 


30 


/35 


40 


PROGRAM S AMO AT ( INPUT. OUTPUT# TAPE5-IMPUT.TAPE6-0UTPUT) 

C 

DIMENSION A(16)#B (6)»0(16).R(4).ATIL(16)#BTlL(6).QTIL(16).Wm(8). 
1RTIL(4).FTU(8).AHAT(16).QHATU6).PHAT(8).0UMMY(116)»P( 16).F(8) 
DIMENSION NA(2).N6(2)fNQ(2).NR(2)»NATIL(2).NBTll(2).NQTIL(2)»NWTU 
1(2).NRTIL(2).NFTIL(2).NAHAT(2).NQHAT(2).NFHAT(2).I0P(5).NP<2).NF(2 
2).N0UM(2) 

C 

LOGICAL lOENT.OISC.NEWT. STABLE* FNULL 
C 

c 

C INPUT HOLLERITH DATA FOR TITLE OF OUTPUT 
CALL RDTITL 
C 

C INPUT COEFFICIENT AND WEIGHTING MATRICES FOR CONTINUOUS SYSTEM 
CALL READ(4*A*NA* 6*NB*0*N0*R*NR) 

C 

C 

C GENERATE EXP (A«0ELT) AND (INTEGRAL EXP ( A^TAU) ) ^'B 
0ELT-.05 
I0P(1)«0 

N1 « (NA(1)A*2)^1 
C 

CALL EXPINTCA.NA* ATIL.NATIL.OUMMY.NDUM.DELT* lOP* DUMMY (Nl) ) 

CALL MULT(0UMNY*N0UM*S*NB*BTIL*N9TIL) 

CALL PRNT(ATIL*NATIL*4HATIL.1) 

CALL PRNT(BTIL*NBTIL*4HBTIL*1) 

C 

C 

C GENERATE DIGITAL PERFORMANCE INDEX WEIGHTING MATRICES 
10P(2)-1 

CALL EQUATE(0*N0, OTIL.NQTIU 
CALL £OUATc(R*NP,RTIL*NRTIU 

CALL SAMPL(A*NA*B.NB.0TIL.NQTIL*RTIL*NRTIL.HTXL.NWTIL.0ELT. IDP.OUH 
IMY) 

CALL PPNT(0TIL*NQTIL*4HQTIL*1) 

CALL PRNT(WTIL.NWTIL#4HWTIL.1) 

CALL PRNT(RTIL*NRTIL*4HRTIL#1) 

C 

C 

C FIND PREFILTER GAIN WHICH ELEMINATES CROSS-PRODUCT TERM 
C IN DIGITAL PERFORMANCE INDEX 


U1 

U1 



PP06RAM SAHOAT 


7<»/74 OPT-1 


FTN 


77/07/27* OS«36«31 


45 


50 


55 


60 


65 


70 


75 


10P(3)»1 

CALL EQUATE (ATXL»NAril»AHAT»NAHATI 
CALL EQUATE (QTIL#NQTlt#QHAT# NOHAT) 

CALL PPEFlL(AHAT»NAHATf BTlL>NBTIt#OHAT#HQHAT#WTTLfNWTIL#RTIL#HRTIL 
1»FTIL#NFTIL> IOP*OUH^Y) 

CALL PPNT(AHAT#NAHAT»4HAHATf 1) 

CALL PRNT(0HAT#N0HATf4HQHAT>l) 

CALL PRNT(FTIL#NCTILf4HFTIL#l) 

C 

C 

C SOLVE FOP F HAT 
IDENT- .TRUE. 

DISC • .TRUE. 

NEWT ■ .TRUE. 

STABLE • .FALSE. 

FNULL ■ .TRUE. 

ALPHA • .9 
lOPd) • 1 
I0P(2) • 0 
ICIP(3) - 0 
I0P(4) • 1 
I0P(5) • 1 

CALL ASYPREG(AHAT»NAHAT. BTIL»NBTIL^H»NHf QHATfNOHAT.RTUfNRTIL^FHAT 
I»NFHAT>P»NP. IDENTf DISC# NEWT. STABLE# FNULL# ALPHA. XQP.DUHHY) 

CALL ADD(FTIL#NFTU#FHAT#NFHAT#F#NF) 

C 

C 

C OUTPUT F AND F HAT GAINS 
CALL LNCNT(IOO) 

PRINT 100 

100 F0RMAT(9 for THE ORIGINAL SAHPLED^DATA PROBLEHA) 

CALL PRNT<FHAT#NFHAT#4HFHAT#1) 

CALL PRNT(F#NF#4HF #1) 

C 

C 

STOP 

END 



EXAMPL6 2 - OPTIHAL SAMPLED^OATA REGULATOR 


A MATRIX 

-2*6000000r:+00 

-7.5000000E-02 - 

7.8000000E-02 
l.OOOOOOOEtOO 

8 MATRIX 
l»7000000t+01 
8.2000000E-01 

0 . 

0« 

Q MATRIX 
l.OOOOOOOE+00 

0. 

0. 

0. 

R MATRIX 

1,0000000£+00 

0 . 

ATIL MATRIX 
0.746O216E-O1 
-3.06^9872£-03 
3.77A3137E-03 
A.602O759E-O2 

BTIL MATRIX 
7.9692392E-01 
3,925021PE-02 
: 6,180065AE-0«» 
2»0^35169E-02 

OTU MATRIX 
A,3959<»0AE-02 
e.73987A0E-C<r 
-4.1390797E-02 
1.160879AE-03 


G3 


<» ROWS ** COLUMNS 

2.50000006-01 -3.8O0UO0OE+OI 
2*7^000006-01 <».AOOOOOOE^OO 
9.i9QOOOOCE-01 -2. 3000000E-01 
7.3000000F-02 0. 

V ROWS - 2 COLUMNS 

7.0000000 E+oa 
3.2000000E+00 
A.6000000E-02 

0 . 

A ROWS 4 COLUMNS 

0 . 0 . 

l.-O'OOOOOO.E+OO 0. - 

0. i.'Oooooooe+oQ 

0 . ‘ , 0 . 

2 ROWS 2 COLUMNS 

Q. 

l.OOOXtOOOF + 00 

• ROWS ^ COLUMNS 

5.62065>53'6.-02 -1 . 76<i A03 5E +00 

9.811^6906-01 2.1998505E-01 

A. 87079236-02 9. 795 8A3 7E-01 

A.9177112E-03 -A. A 80 976 8 E-02 

A ROWS 2 COLUMNS 

3.223A606E-01 
1. 58957336-01 
6..ab76l.37F-03 
7.98528A5E-03 

% Raws A COLUMNS 

«.73987,AOE^O^. --A. 1 39079 7^-02 
A .9 ZZ 7 561-E - 02 . 2.861 932 7 E-D3 

2.8619327€rf03> . 1 *0-38 67*23 £ -01 
1.1^66907E-0A -6. 383 598 1 E-OA 

A ROWS 


WTIL MATRIX 


2 COLUMNS 


ORACLS PROGRAM 


0. 

0 . 

5.2000000E-02 

0 . 


0. 

0. ‘ 

0 .' 

l.OOOOOOOE+00 . 


-2.352A289E-03 

2.8616A68E-0A 

2.5772852E-03 

9.9996069E-01 


l.l608s79AE-03 

1.1266907E-0A 

-6.3835981E-0A 

A.9999173E-02 



Ul 

00 


EX4HPLF c - nPTIP4L SA^fPlED-OATA REGULATOR ORACLS PROGRAM 

3.7319029E-02 1 , 5229661E-02 

3.1985357E-03 -7. 38388 7P E-03 

-4.8263666E-02 -2.06678 55E-02 

6.A086387E-0't 2.^9911^'.E-0‘« 

RTIL MATRIX 2 ROWS 2 COLUMNS 

6.096'(766E-02 A . 39 519 93 E-03 
A.3A51993E-03 5. 22 3A106E-02’ 

ahat matrix a rows a columns 

5.985A577E-01 5.A76A99AE-02 -1.A0A9026E+00 -7.0633127E-03 

A. A95A715E-03 9 . 68 26103 E-01 2.0868737E-01 3 . 97980A lE-OA 

2.7592A06E-03 -A . 82 2398AE-02 9. 8096098E-01 2 . 560617AE-03 

3.9775917E-02 A. 8601B65E-03 -3. 5637777E-02 9.998A0A3E-01 

QHAT matrix a rows A COLUMNS 

3.7A87559E-G2 8 . A5 165 8 0 f -OA -3 .2 96212 6E-02 1 .050A AA 5E-03 

8.A516580E-0A A . 39066 A2 E-02 2. 8 63 593 2E-03 1 . 117A823E-0A 

-3.2962128E-02 2. 8635932 t-03 9.238606AE-02 -A. 9A58179E-0A 

1.050AAA5E-03 1.117A823E-0A -A. 9 A56 179E-0A A . 9997293E-02 

FTIL matrix 2 ROWS A COLUMNS 

2.97AA352E-01 3 , 1 A5686 1 E-02 -3. 8A00852E-01 5 . 1 1 58A7AE-03 

1.2103929E-C1 -7 .32 975 16E-02 -1. 656 9A1 8 E-01 1.9666527E-03 


COMPUTATION OF F SUCH THAT 4-8F IS ASYMPTOTICALLY STABLE IN THE DISCRETE SENSE 


A HATRIX A ROWS A COLUMNS 

6.6505086C-01 6.08A9993E-02 -1. 5610029E+00 -7. 6 A812 52E-03 

A.99A9683E-03 1 .0758A 56 t +00 2. 3 187A8 6E-01 A . A2200A6E-0A 

3.0653229E-03- -5 . 3 5 32205 t-02 1.0899566E+00 2 . 6A5130 AE-03 

A.A195A6AE-02 5. A002072E-03 -3. 9 597530E-02 1 . 1109338E+00 

8 MATRIX A ROWS 2 COLUMNS 

8.85A7102E-01 3 . 5 81 62 2 9t -01 

A.3611353E-02 -1 . 766192 5 E-01 

6.3667393E-0A 7. 63068 19E-03 

2.27057A3F-02 8 . 37253 »3F-03 



EXAMPLE 2 - OPTICAL SAMPLEO-OATA REGULATOR ORACLS PROGRAM 


ALPHA • .33873839E+00 


F MATRIX 2 ROWS A COLUMNS 

8.1A98121E-01 2*2921930E + 00 •3.0628590E + 01 1. A89992AE4'01 

1.0955671E-01 -5. 7A22787E +00 8, 0373720E+01 -3, A265557E+00 

A-BF MATRIX A ROWS A COLUMNS 

-9.5830A73E-02 9.78A72 50E-02 *3. 2271096E+00 -1 . 197A036E+01 

-l#11976AU-02 -3.9317030F-02 1.5763176F+01 -1.25A5593E+00 

l*67020A2E-03 -1 . 1 33 86 92 E-02 A, 976622 lE-01 !• 8760698E-02 

2.A7l8663t-02 A . 3028601 E-03 -5. 7271536E-02 8 . 0302 221E-01 


EIGENVALUES OF A-BF 

EIGN MATRIX A ROWS 2 COLUMNS 

2. 3A 192296-01 3, 35 28007E-01 

2, 3A 192296-01 -3. 3 5280C7E-01 

3*A908617E-01 3 . 1 66A5 98 E-01 

3.A908617E-C1 -3.166A598E-01 


MODULI OF EIGENVALUES OF A-BF 


MOO MATRIX A ROWS 1 COLUMNS 

A. 08972816-01 
A.0897281E-01 
A,7130228c-01 
A.713022»F-01 


PROGRAM TO SOLVE DISCRETE STEADY-STATE RICCATT EQUATION 8Y THE NEWTON ALGORITHM 


A MATRIX pows A COLUMNS 

5.985A577E-01 5 . A76A99AE-02 -!• A0A9026E+00 -7,C633127E-03 

A.A95A715E-03 9.68261036-01 2.0668737E-01 3 .97980A1 E-OA 

2.7592A06E-03 -A. 822398 AE-02 9. 8096098E-01 2. 560617AE-03 


U1 

VO 


CTn 

O 


EXAiiPLE 2 - OPTIMAL SAMPLEO-OATA REGULATOR 


0RACL5 


3.Q775917E-02 


0 MATRIX 
3.7487559E-02 
8.A516580E-04 
-3*2962l28E-02 
1.050^4A5E-03 


A*8601665E-03 


^ ROWS 

8*45165806-04 

4.8906642E-02 

2.8635932E-03 

1.1174823E-04 


-3. 56377776-02 
2 COLUMNS 


4 COLUMNS 
-3.2962128E-02 
2.8635932E-03 
9.2986064E-02 
-4,9458l79£-04 


9.9984043E-01 


1*05044456-03 

1*11748236-04 

-4.94581796-04 

4.9997293E-02 


8 MATRIX 4 ROWS 

7.9692392E-01 3.2234606E-01 

3*92502186-02 -1*58957336-01 

6.1800654E-04 6* 8676137F.-03 

2.0435169E-C2 7* 9852P 456-03 


H IS AN IDENTITY MATRIX 


R MATRIX 2 ROWS 2 COLUMNS 

6.0964766E-02 4.34519936-03 

4 *34519936-03 5*22341066-02 

INITIAL F MATRIX 


F MATRIX 2 ROWS 4 COLUMNS 

0* 14981216-01 2*29219306+00 -3*06285906+01 

1 *0955671 E-01 -5 . 7422787E +00 8*03737206+01 


1.4899924E+01 

-3.4265557E+00 


10 ITERATIONS TO CONVERGE 


FINAL values OF P AND F AFTER 


P MATRIX 4 ROWS 

5.1520273E-02 1*05909666-02 

l*0590966E-02 3.68958166-01 

-9*50897216-02 -3.4132485E-01 

5 *61737306-02 4.2019O65E-O2 


F MATRIX 
2.7657460E-01 
8*24146746-02 


2 ROWS 

4*45897476-01 

-9*45968766-01 


4 COLUMNS 
-9*50097216-02 
-3*41324856-01 
2.92110596+00 
-1*34984716-01 

4 COLUMNS 
-1*56785306+00 
4*63133306-01 


5*61737306-02 

4.28190856-02 

-1*34984716-01 

1.0635575E+00 


6*75154376-01 

1.38775416-01 


PROGRAM 



EXA^PLP 2 - 0*>TIWAL SAHPLEO-DATA REGULATOR 


RESIDUAL ERROR IN RICCATI EQUATION 


EROR NATRIX A ROWS A COLUMNS 

-3«10862456-15 -4.2188^»75E-15 1. 1990A09E-1A 

-A*2743586E-15 -7, 105427AE-15 !• 776356SE-1A 

l«1990A09t-lA 1.7763568E-1^ -4. 2 63 2564E-1A 

-A.6851A12E-1A -7. 5051076fc-l^t 1. 9895197E-13 


EIGENVALUES 0^ P 


EVLP MATRIX ROWS 1 COLUMNS 

4.5717003E-02 
3.2330141E-01 
1.0566117E+00 
2V97951TTE>00 


CLOSED-LOOP RESPONSE MATRIX A-8F 


•A-3F MATRIX A ROWS A COLUMNS 

3. 515708 lE-01 A.3A79A0AE-03 -3. 0A7322 AE-01 

6.7A027AAE-03 8. 00390 79E-01 3. A38 AA39E-01 

2*0223236E-03 -4.200300A E-02 9. 787493 IE-01 

3.3465964E-02 3.3020260E-03 -7.2966874E-03 


eigenvalues of A-BF 


3.8573552E-01 0. 

8#5933093E-01 7.9976165E-02 

848933093C-01 -7.99761 85E-02 

9.5124891E-01 0. 


ORACLS PROGRAM 


-4.6185278E-14 

-7.5495166E-14 

1.9895197E-13 

-8.3133500E-13 


-5.6984369E-01 

-4.0426071E-03 

1.1903116E-03 

9.8493537E-01 



O' 

ro 


EXAMPLE 2 - optimal SAMPLED-OATA REGULATOR ORACLS 

FOR THE ORIGINAL SAMPLED-OATA PROBLEM 


FHAT MATRIX 2 ROWS 

2.7657A60E-01 A,A5S97A7E-01 

8.2A1A67AE-02 -9.A596876E-0I 


F MATRI 

5.7A01812E-01 
2.03A5396E-01 


X 2 ROWS 

A.7735A33E-01 
-1.01R2663E+00 


A COLUMNS 
-1.5678530E+00 
A.6313338E-01 

A COLUMNS 
-1.9516615E+00 
2.9723921E-01 


6.7515A37E-01 

1.38775A1E-01 


6.8027022E-01 

1.A07A207E-01 


PROGRAM 



Example 3 - Model Following 

This problem demonstrates the model-following capability of ORACLE and 
illustrates the output format of the transient response subroutine TRANSIT. 

Given the linear system 

x(t) = A x(t) + B u(t) 


with x(0) = xq given and A and B as stated in example 1, use the explicit 
model-following procedure to find a control law which causes the state x(t) 
to track the state Xnj(t) of the system, 


~ Ajjj Xjj|(t) + Bjjj Ujjj(t) 

where Xn,(0) = Xj„ is given and the step input is 


UfflCt) = 


2 

0 


with 


-0.981 

0.177 

-10.0 

0.0 

0.030 

-0.092 

5.23 

0.0 

0.0 

-1 .0 

-0.732 

0.052 

1 .0 

0.0 

0.0 

0.0 


6.3^ 

4.58 

0.690 

-2.65 

0.0 

0.0 

0.0 

0.0 


Generally, the state and model dynamics represent two sets of linearized equa- 
tions and u(t) is to cause the system represented by x(t) to perform like 
the model Xnj(t) when the model is subjected to a particular step input. Data 
for the model dynamics are from handling-quality model 3 of reference 31. 
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The solution to the posed problem can be found by using the subroutine 
EXPMDFL of ORACLS. The foregoing problem can be placed in the format of 
EXPMDFL by redefining the model dynamics as 


*m^^^ “ ^m *m^^^ 



and 

where 




where I 4 is the x 4 identity matrix. Weighting matrices in EXPMDFL are 
chosen to be 

Q = diag ( 10 . ,10. ,10. ,10.,) 
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and R is the 2x2 identity matrix. A transient response is evaluated 
using TRANSIT with 


x(0) = x^ = 0 


The executive program and output data follow. 
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ON 


PROGRAM MOOrOL OPT-1 


FTN 4*6+439 77/07/27. 08«36#41 


1 


5 


10 


15 


20 


25 


30 


35 


40 


PROGRAM MOO FOK INPUT# OUTPUT# TAPE 5» I NPUT/TAPE6*0UTPUT) 

C 

DIMENSION A(16)#B(8)#AM(36)#HM(24)#Q(16)#R(4)#F(20)#P(40)#0UMMY(70 
10)#AC(100)#9C(20)#X(10) 

DIMENSION NA(2)#NB(2)#NAM(2)#NHM(2)#NQ(2)#NR(2)»NF(2)#NP(2)#NAC(2) 
1>NBC(2)#I0P(5) #N0UM(2J#T(2)#NXC2) 

LOGICAL DI SC# NEWT# STAB LE#HI0ENT#HMID6NT 
C 
C 

C INPUT HOLLERITH DATA FOR TITLE OF OUTPUT 

CALL RDTITL 
C 

C INPUT COEFFICIENT MATRICES FOR PLANT ANO MODEL 

CALL READ (5# A# N A# B# NB# AM# N AM# HM# NHM# X# NX ) 

C 

C define weighting MATRICES 

N0(1)» NA(1) 

N0(2)* NO(l) 

NP(n« NB(2) 

NR(2)« NR(1) 

CALL UN1TY(0#NQ) 

CALL SC'ALE(0#NQ#Q,NQ#10*) ^ ... ... . 

CALL UNrTY(RfNR) 

C 

C 

C COMPUTE MODEL-FOLLOWING GAINS 

lOP(l) • 1 
I0P(2) • 1 
l0P(3)-0 
I0P(4) • 0 
ICP(5) « 1 
DISC • .FALSE. 

NEWT • .FALSE. 

STABLE • .false# 

HIOENT-.TRUE. 

HMTDENT ■ .FALSE. 

C 

CALL EXPMDFL(A#NA#B#NB#H#NH#AH#NAM#HM#NHM#Q#NQ#R#NR#f#NF# P#NP#H10E 
1 NT# HMICFNT# DISC# NEWT# stable# FNULL#ALPHA#IOP# DUMMY) 

C 

c 

C GENERATE COEFFICIENT MATRICES FOR INPUT TO TRANSIENT 



PROGRAM MOOPOl 


7<»/7<» DPT«1 


FTN 4«6^439 


77/07/27. 08.36.A1 


A5 


$0 


55 


60 


65 


70 


75 


C RtSPOMSE SUeROUTIKE 
NDUM(l) • NAM(l) 

N0UM(2) • NAU) 

CALL NULL(DUMMY.NDUM) 

CALL JUXTR(A#NA#DUMMY.NOUM#ACfNAC) 

NDUM(l) • NA(1) 

NrUM(E) . NAM(1> 

CALL MULLCOUMMY.NOUM) 

N1 • NDUM(1)+NDUM(2) ♦ 1 

CALL JUXTR( DUMMY. NOUM. AH. NAM# DUMMY(N1)*NHM) 

N2«NHM(1)*NHM(2)+N1 

CALL JUXTC(ACf NACf DUMMY (MU# NHM. DUMMY (N2)#N0UM) 

CALL SOUATEtDUMMY (N2 ) #NOUM# AC #NAC I 
NOUM(l) ■ NAM(l) 

N0UH(2) - NB(2) 

CALL NULl (DUMMY. NOUM) 

CALL JUXTR(B#NB. DUMMY# NDUM#BC#NBC) 

C 

C COMPUTE TRANSIENT RESPONSE TO OBSERVE MODEL-FOLLOWING ACCURACY 
IQP(l) - 0 
I0P(2) ■ 0 
I0P(3) • 0 
T(l) • 50. 

T(2) • 2. 

C 

CALL LNCNT(IOO) 

PRINT 100#NA(1).NAM(X) 

100 FORMAT!* IN THE TRAJECTORY OUTPUT TO FOLLOW THE FIRST*IA* CO 

ILUMNS CORRESPOND TO X TRANSPOSE*/* AND THE NEXT*U* COLUMNS TO X 
2M TRANSPOSE*) 

C 

CALL TRANSIT ( AC. NAC. BC. NBC# H.NH.G.NG.F.NF#V.NV#T.X# NX# DISC# STABLE# 
1I0P#DUMMY) 

C 

C 

STOP 

END 





; O IMOOOOOOOO O O O M OOH'OUi-i) 


6XA*1PU6 3 - MOOEL-FCLLOWING 


ORACLS PR0GRAH 


A MATRI] 

-2#6000000t+00 
-7.5000000E-02 
7.8000000E-02 
l.OOOOOOOE+00 

B MATRI 

1.7000000E+01 

8,2000000e-01 

0. 

0« 

AH HATRI 
01OOOOOE-O1 
0000000 E-02 

OOOOOOOE+00 


A ROWS 

2.5000000E-01 

-2.7000000E-01 

-9.9000000E-01 

7.80p00p0E-02 

: A ROWS 

7.0000000E+00 
-3.2000000E+00 
A,6000000E-02 

0, 

( 6 ROWS 

1.770000CE-01 
-9.2000000E-02 
-l.OOOOOOOE+00 

0. 

0* 

0. 


A COLUMi^S 
-3*8000000E+01 
A.A.OOOOOOE+00 
-2.300000.0E-01 

0 . 

2 COLUMNS. 


6 COLUMNS ■ 
-l*OOOOOOOE>01 
5#2300000E+00 
-7.3200000E-01 

0 . 

0* 

0* 


0« 

0 *.. 

5.2JOOOOOOE-02 

0 . . 


0* 

0 • 

5*2000000E*02 

0. 

0. 

0* 


6«3AOOPOj0Et,OO 
-6.9000000E-01 ^ 

0. 

0. 

0« 

0. 


A.saoooooE^oo . 

-2.6500000E+00 

0 . 

0. 

0* 

0 . 


HM MATRIX 

A ROWS 

6 COLUMNS 


oooopooe+ 0,0 

0« 

0 . 

0* 


UOOOOOOOE+00 

0. 

0* 


0* 

l.OOOOOOOEAOO 

0* 


0* . 

0. 

l.OOOOOOOE^OO 

X MATRIX 

10 ROWS 

1 COLUMNS 



•OOOOOOOE^OO 



EXANPie 3 • iOOEl-FOUOtfISG 


ORKLS PROGRAN 


fROGRAH TO SOLVE ASYMPTOTIC CONTINUOUS EXPLICIT NOOEL-FOLIOWING PROSLER 
PLANT OYMANICS 


A MATRIX 

<t ROWS 

4 COLUMNS 



:-2*6000000E+00 

2.5000000E-01 

-3.8000000E+01 

0. 


-7#5000000E-02 

-2.7000000E-01 

4.4000000E«00 

0. 


7t8000000£-02 

-9*9000000fc-01 

-2.3000000E-01 

5.2000000E-02 


l.OOOOOOOE^OO 

7.80000006-02 

0. 

0. 


B MATRIX 

9 ROWS 

2 COLUMNS 



1*7000000E>0I 

7.0000000E+03 




$•20000006-01 

-3.200000C6+00 




0« 

<»»60000006-02 




0« 

0« 




H JS AN IPENTITY. 

MATRIX 

. , 

. . 

- •, 


NOOEl DYNAMICS 


AM Matrix 
-9#$ioooooE-oi 

6 ROWS 

1*77000006-01 

6 COtUMMS 
-l.OOOOOOOE+Ol 

0* 

6*34000006400 

4*58000006400 

a.OOOOOOOE-02 

-9*20000006-02 

5*23000006400 

0* 

-6*90000006-01 

-2*65000006400 

0. 

-1. 00000006+00 

-7.32000006-01 

5*20000006-02 

0* 

0* 

1*00000006400 

0. 

0* 

0* 

0* 

0* 

0% 

0. 

0* 

0* 

0* 

0* 

0« 

0. 

0* 

0* 

0« 

0« 

HM MATRIX 

X.0000000E400 

4 ROWS 

0« 

6 COLUMNS 

0* 

0* 

0* 

0* 

0» 

l.OOOOOOOE+OO 

0* 

0* 

0* 

0* 

0« 

0. 

1*00000006400 

0* 

0* 

0* 


0* 

0* 

1*00000006400 

0* 

0* 



o 


EXAHPCE 3 - ^fODEL-FOUOWIMG 


ORACLS 


weighting matrices 


Q MATRIX ^ ROWS 

l.OOOOOOOE+01 0. 

4 COLUMNS 

0. 

0. 

0. 

l.OOOOOOOE+Ol 

0 1 

0* 

0. 

0. 

l.OOOOOOOE+01 

0. 

0. 

0. 

0. 

l.OOOOOOOE+01 

R MATRIX 2 ROWS 

1.00000002400 0. 

0. 1.0000000E400 

CONTROL LAW U • -Ft COL.tX.XH) ) 

PART OF F multiplying X 

Fll MATRIX 2 ROWS 

2.8715216E+C0 1 . 1 8677 8R E +00 

2 COLUMNS 
. F . (F11,F12) 

A COLUMNS 
-2.3A93987E+00 

3.0218811E+00 

1,0990931E+C0 

*3.03621^1E+0Q 

1.592190AE+00 

9.3351305E-01 

Pll MATRIX A ROWS 

1.678A659E-01 2 . 21091 23E-02 

A COLUMNS 
-1.1051966E-01 

1.7350075E-01 

2.2109123E-02 

9.80O31«»OE-O1 

-5. 73859108-01 

8.825^03<rE-02 

-1.1051966E-01 

-5.7385910E-01 

1.1510A12E+01 

3.088^697E-02 

1.7350075E-01 

e.325^03AE-02 

3.088A697E-02 

1.0164^7<»E+01 


eigenvalues of pu 


X,63^60<»6E-O1 
9.57109<»26-01 
1.01679<»2E + 01 

PLANT CLOSEO-IOOP RESPONSE MATRIX A - 8F11 


RROGRAN 



example 3 


MOOEt-FOLLOWiNG 


ORACLS PROGRAM 


-5.9109518£+01 

1.0874501E+00 

2#7441719E-02 

l«OOOOOOOE^OO 


1«3282587E400 

-1.0959044E^0l 

-0.5O33415E-O1 

7.0OOOOOOE-O2 


-9,2055556E^00 

1.1421516E+01 

•3%0324076E-01 

0 . 


-5.7906570E+01 

5.0929927E-01 

9.0583996E-03 

0 . 


EIGIMVALUES OF CLOSED-LOOP RESPONSE MATRIX 


-1.0076624E+00 0. 
-1.2943964E+00 0* 
-9.9317054E+OO 0. 
-5.8137959E+01 0. 

PART OF F MULTIPLYING XM 


F12 MATRIX 
-2#9990036E+00 
-9.50610406-01 


2 ROWS 

1.0535303E^00 

2.9377885E400 


6 COLUMNS 
9,7142903E-01 
*1.5695172E^00 


-3.0178671E+00 

-9.2056907E-01 


-5.5004933E-01 

8.1556250E-01 


P12 MATRI 
-1.7290904E-01 
-7.2521045E-O2 
5.9484706E-01 
-1.5393360E-01 


4 ROWS 

-1.6O50975E-O2 

-9.5186308E-01 

9.2161095E-O2 

1.4943541E-01 


6 COLUMNS 
3.69716736-02 
4.18103646-01 
-1. 06550266+01 
4.0497750E-01 


-1.7312936E-01 

-9.1058595E-02 

-1.1101935E-03 

-1.0136463E+01 


-2.0297024E-02 

-2.4999989E-01 

3.4270003E+00 

6.3234900E-02 


-4.3236218E-01 

4.4887921E-01 


-1.9460473E-02 

-1.2382211E-01 

4.1059083E+00 

-3.7895492E-01 



EXAMPLE 3 • MaOEL-FOLLOWlNG 


ORACLS PROGRAH' 


IN THE TRAJECT3RY OUTPUT TO FOLLOW THE FIRST ^ COLUMNS CORRESPOND TO X TRANSPOSE 
AND THE NEXT 6 COLUMNS TO XH TRANSPOSE 


COMPUTATION OF TRANSIENT RESPONSE FOR THE CONTINUOUS SYSTEM 


A MATRIX 

2.6000000E+00 

10 ROWS 
2.5000000E-01 

10 COLUMNS 
-3.8000000E+01 

0. 

0. 

0. 

0« 

0. 

7.5000000E-02 

0. 

-2.7000000E-01 

0* 

A.AOOOOOOE+00 

0. 

0. 

0. 

0. 

0 • 

7.8000000E-02 

0. 

-9.9000000E-01 

0. 

-2.3COOOOOE-01 

5.2000000E-02 

0« 

0* 

0. 

0* 

l.OOOOOOOE+00 

0 * 

7.S000000E-02 

0. 

0* 

0. 

0. 

0. 

0. 

0. 

0. 

0* 

0. 

0. 

0. 

0. 

-9.8100000E-01 

1.7700000E-01 

-I.0000000E401 

• 

o 

• 

o 

6.3‘t00000E«00 4.5B00000E4’0C 

0. 0. 

1 

0. 

3«0000000E-02 

-9.2G00000E-02 

5.2300000E400 

0. 

0. 

-6.9000000E' 

0. 

-01 -2#6500000E+OC 

0. 

1 

0. 

0 . 

-1«0000000E400 

-7.3200000E-01 

5.2000000E-02 0. 

0. 0« 

0. 

0. 

0. 

1.0000000E400 

0. 

0 . 

0. 

0« 

• 

o 

• 

o 

• 

o 

• 

o 

0. 

0. 

0. 

0. 

• 

o 

• 

o 

o 

• 

o 

• 

0. 

0« 

0* 

0» 

0» 

0. 

0. 

0. 

0. 






B MATRIX 10 ROWS 2 COLUMNS 

1.7000000E401 7.0000000E+00 

8.2000000E-01 -3.2000000E+00 

0. A.6000000E-02 

0 * 0 » 

0 « 0 * 

0 * 0 » 

0* 0. 

0« 0« 

0 . 0 . 

0* 0* 


H IS A NULL MATRIX 



£X*»(PLE 3 - MODEl-FOtlOWING 


ORACLS 


PROGRAM 


6 IS A NULL MATRIX 

F MATRIX 2 ROWS 10 COLUMNS 

2.8715216E+00 1. 1967788E+00 -2.3A93987E+00 3.0218811EA00 -2.9990036E+00 -1.0535303E+00 9.71A2903E-01 

-3.0178671E+00 -5.5009933E-01 -A. 3236218E-01 

X.0990931E+00 -3.0362141E+00 1.592190AE+00 9.3351305E-01 -9. 50610AOE-01 2.9377885E+00 -1.569S172E«00 

-9.2056907E-01 8.1556250E-01 9.98879212-01 


V IS A NULL MATRIX 


COMPUTATION OF THE MATRIX EXPONENTIAL £XP(A T> BY THE SERIES METHOD 


A MATRIX 10 ROWS 10 COLUMNS 


■5.9109518E9CI l.3292587t VOO -9.2055556E+00 -5.7906570E+01 

5.7797725E+01 3.6919010S+00 9.2080025E+00 

5.7637333€m 

-2.65950^%E^00 

-5«3276735E900 

1.0879501E+00 -1.0959099E+01 1.1921516E+01 5.0929927E-01 

-9.7116998E-01 3.0608909E+00 1 .7909505E+00 

-5.8277036E-01 

l«026^818Em 

•5*619O267E^00 

2.7991719E-02 -8.5033915E-01 -3.0329076E-01 9.0583996E-03 

9.2396177E-02 -3.751i675E-02 -2.0698999E-02 

4.372807SE-02 

-1#3513827E-01 

7*2l97789E-!02 

1*OOOOOOOE900 7.8000000E-02 0« 0* 

0* 0. 0. 

0. 

0. 

0 . 

0. 0. 0. 0. 

0* 6.3900000E+00 9. S800000E900 

-9.8100000E-01 

1«7700000E-01 

-ItOOOOOOOE^Ol 

0. 0. 0. 0# 

0. -6.9000000E-01 -2. 6500000E900 

3«OOOOOOOE-02 

-9«2000000E-02 

9.2300000EV00 

0. 0. 0. 0* 

$•20000002-02 0* 0. 

0 . 

-1«0000000E^00 

-7*3200000EH>1 

• 

O 

• 

o 

• 

o 

• 

o 

• 

o 

• 

o 

• 

o 

l.OOOOOOOE^OO 

0. 

0« 

o 

• 

o 

• 

o 

• 

o 

• 

o 

• 

o 

• 

o 

• 

0. 

0. 

0* 

• 

o 

• 

o 

• 

o 

• 

o 

• 

o 

• 

o 

• 

o 

0. 

0 . 

0 % 


T • 


20000C00E901 


example 3- model-following 

EXPA MATRIX 10 ROWS 10 COLUMNS 


ORACLS PR06RAH 

■2*38911316-03 -!• 

1»2605585E-01 

3251850E-03 -1. 138 5082E-04 -1#3738924E-01 

3.9547601E+00 -1*4409641E400 

1.19556056-01 

1.06253376-01 

2.14275256+00 

1.3074054E-0A -8, 

5.3247219E-02 

56705 45E-03 9. 7002193E-02 -8 • 5196158E-04 

7.2397266E-01 3. 8807667E-01 

6.26447296-02 

1.45347876-02 

-1.05993496+00 

1.4763163E-05 -7. 

-l«6411698E-02 

2188507E-03 8*2534175E-02 6# 5207081E-03 

8 •236484 5E-02 4* 083 1924E-01 

9.21765346-03 

2.22590466^01 

-1.55152216-01 

2.3921962E-03 2* 

7.6035997E-01 

29U988E-03 -1.089844 0E-02 1 • 3680976E-01 

5.85625506+00 2 . 3683158E-01 

8.48213866-01 

2.00373916+00 

-1.22955796-01 

0* 

-7.8748841E-03 

0. 0. 

4.O115443E+0C -1.44220676+00 

1.21072986-01 

1.16708756-01 

2.13815116+00 

0* 0« 

5.4105409E-02 

0. 0. 

3.07012146-01 1.17170376-01 

4.80069116-02 

-3.54541426-02 

-1.04450386+00 

0« 0* 

-9.3801938E-O3 

0. 0 . 

1.77194266-01 5.77030786-01 

1.28708456-03 

2.00669176-01 

-1.16933786-02 

0* 0. 

0.9663838E-OI 

0. 0. 

5.09262936+00 2.63047126-01 

8.52275166-01 

2.01721776+00 

-1.51440086-01 

• 

o 

• 

o 

• 

o 

0. 0. 

1.00000006+00 0. 

0. 

0. 

0. 

• 

o 

• 

o 

• 

o 

0. 0. 

0. 1.00000006+00 

0. 

0. 

0. 


structure of printing rn follow 

TIME OR STAGE 

state - X TRANSPOSE - FROM OX » AX ♦ BU 

OUTPUT - Y TRANSPOSE - FROM Y « HX ♦ GU IF DIFFERENT FROM X 
CONTROL - U TRANSPOSE - FROM U « -FX ♦ V 


0 . 

0 « 0 * 0 « 0 » 0 # 0 * 0 « 

2.00000006400 0. 


0 « 



EXAMPLE 3 - rtOOEL-FOUOWING 


1.10009S7E4C0 f-1.6311250E^00 

•ZOOOOOOE^Ol 

7.9095202E^00 1.4479453E^00 !• 6472969E-01 

1.1785259E+01 2 .OOOOOOOE^OO 0. 

l,5927996Ei‘00 1.0A10762E-01 

•AOOOOOOH^Ol 

9.54893X7E+00 2.2039396E+00 2.0659338E-01 

3.0375205E^0l 2*OOOOOOOE+O0 0* 

1.8227l^8E^00 1.0305996E-01 

.6000000E+01 

9#6966669E^00 3.271631<^E400 1.7391236E-01 

A.9769107E+01 2.0000000E^OO 0. 

1.79A2873E+00 6.0926728E-02 


--4 

VJ1 


ORACLS PR06RAN 


l«l712510E<t>0l 8»0230886E400 6«1^02928£«*01 3«5^3888lE-01 


3.0288443E+01 9. 7310589E+00 1.29A90A9E^00 3«7319936E-01 


9,9680598E*01 9.9051858E+00 2.2907575E+00 3.2n9551E-01 



-<3 

o\ 


EXAMPLE 3 - 

#8000000E4>01 

9«A569500E>00 < 

6 « 94234996>01 

1«8504096E400 3 

•1000000E402 

9«397<i701E400 ! 

8.9105733E+01 

l.90971^2E+00 

•1200000E^02 

9.3^39028E+00 i 
l#0876250E+02 

1#9994176E400 - 


MOOEL-POILOWING 


ORACLS PROGRAM 


•3502637E400 2* 1320376E*0X 6. 933I980E^01 9«7973579E^00 3.3393A57E900 3«5575370E*01 

2.0000000E+00 Qm 


•1082298E-02 


•3531655E+00 2. ^666911E-01 8.9012688E401 9#8l53021€+00 4*3490632EtOO 3«8519943£-01 

2.0000000E+00 0# 


•6113310E-03 


•3634006E400 2.6356106E-01 l.086689lEf02 9.8410301E4>00 9.3497710E900 3«9a85896EH)l 

2.0000000E+00 


•5910245E-02 


•1400000E^02 


eXAHPLE 3 


MOOFL-FOLLOWING 


9.2595870E+00 7* 3 979771 E +00 2. 9 39^1 9 6E-01 

l,2842<»Alt + 02 2.0D000006 + 00 0. 

!• 9814221 E+00 -9. 9857506E-02 

•1600000E+02 

9.1761171E+00 8.4098289E^O0 3. 071 511 6E-01 

1«4609427F*02 2#0000000c+00 0* 

2#0230646E*00 -1. 4313242 E-Ol 

•1800000Ef02 

9#0997321E+00 9.42742 67E+00 3.2918291E-01 

1.6776481E+02 2.0000000E+00 0. 

2.0634769E+00 -1 .8702727E-01 

•2000000E^C2 

9»02X.7633E+00 1.0447312E + 01 3. 5081025E-01 

1.8743654B402 2 .OOOOOOOE+00 0. 


ORACLS PROGRAM 


1.2832999E+02 9. 8352651E+00 6.3648203E+00 4.1464140E-01 


1.4799891E+02 9.8323699E^00 7.3756775E^00 4»3374640£-01 


1.6766858E^02 9.8353191E^00 8.3841989E400 4.5170904E-01 


1.6733946E^02 9.8370988E^00 9.3941189E^00 4«6920740E*01 


00 


EXAHPte 3 - MODEL-FOLLOWING 
2.10312126>00 -2. 3097232 F-01 

•2200000E^02 

e.9^276^5E+00 !• 1^67616E^01 3. 7275076E-01 

2.0711089E+02 2*0000000t+00 0. 

2#1^31502E400 -2.7A7«692F-01 

• 2 ^ 00000£+02 

8,06<»2<»8OE + OO 1.2<»87737F+01 3.9<^7290 2E-01 

2*2678753b402 2.0000000E400 0, 

2»1832S40E400 -3. 18618 A36-01 

•2600000E402 

8 •785 3 578 E + OO 1 . 3^079 53E401 <»• 166392 <»E-01 

2#^6^6621E402 2 •OOOOOOOE+OO 0. 

2#2233270E400 -3,62A01O<»e-Ol 


2800000E40? 


ORACLS PR06MN 


2«0701295E402 9.8376666E400 1*0A04A5^E4>01 4«86982A6E-01 


2.2666872E402 9# 8387229E400 l«14»1456AE^0l 5.0461197E-01 


2#4636656E402 9# 8399126E+00 1#2424811E+01 5.2257858E-0X 


EXAMPLE 3 - M03EL-P0U0WING ORACLS PROERAH 

8.7073351Et00 1 .4328345E+01 <r.38S5571E>01 2.660A650E+02 9t840979AEP00 1.3A3S217E«01 9.9034652E-01 

2.661A705E+02 2.0000000E+00 0. 

2.2633689E4'00 -A.063A179E-01 

•3000000Et02 

8.6287953EA00 1 . SSABE A3E+01 A.60A8866E-01 2.857286SE+02 9.8420381EA00 1.9445726EP01 5.5813003E-01 

2.8S83006E4'02 2.0C00000£a00 0. 

2t303A352E+00 -A. 50202 fl6 t-01 

.3200000E+02 

8.5502742t+00 1 . 6569A ACE +01 A. 8 2A 2098E-01 3 . 05 A1 296E+02 9.8A31238E+00 1.5A56333E+01 5 . 75913A5E-01 

3.055152AE+02 2.0000000E+00 0. 

2.3A350A3E+00 -A . 9A070 89E-01 
.3A00000E+02 


8.A717A18E+00 1 . 75901 52 E +01 5. 0A3532 5E-01 3. 25099A2E+02 9.8AA2071E+00 1.6A6705AE+01 5.9369631E-01 

3.2520257E+02 2tOOOOOOOE+OC 0. 



CP 

o 


EXAHPL^^ 3 - 
2.3835737E+00 

•3600000E+02 

8i3931956E+00 

3*^«fe9206E+02 

2*^236^86E+00 -I 

•3800000E+02 

8.31^6^21t+00 

3.6^583716+02 

2.<»637285E+00 - 

•4000000E+02 

8.2360809E+00 
3.8<»27752E + 02 

2.5038125E+00 - 


MnOEL-FOLUOWIMG ORACLS PROGRAM 

• 379^^3U-01 


•8610979E+01 5 . 2 62 8 89 5 E-Ol 3 . A97880 5E+02 9, 8952853E+00 1.7977889E+01 6#11^8195E-0X 

2.0000000E+OC 0. 


.81622026-01 


1.9631916E+01 5.A8227A1E-01 3.69<i7883E+02 9.8^63650E+00 1.8^88833E+01 6. 2926992E-01 

2.0000000E+00 0. 


6.2570^506-01 


2.0652965E+01 5. 7016800E-01 3 . 8417 177E+02 9.8474457E+00 1.9499SeBE401 6#4705962E-01 

2.00000006+00 0. 


6.6959191E-01 


4200000E+02 



I 

Example 3 - hooel-following 

8#X575l06E+00 2. 16741 27E ♦©! 5.921 1098E-01 

4.0397349E+02 2.0000000E+00 0. 

2.5439009E+00 -7. 1348413E-01 

• 44000006+02 

8.0789316E+00 2.2695400E+01 6. 1405642E-01 

4*2367162E+02 2 .OOOOOOOE+00 0. 

2.5839937E+00 -7. 57381 14E-01 

•4600000E+02 

8.00034416+00 2. 3716786E+01 6. 3600426E-01 

4.43371926+02 2 .OOOOOOOE+00 0. 

2.6240909E+00 -8 .0126297E-01 

•48000006+02 

7.9217480E+00 2 • 47382 83E+01 6. 5795450E-01 

4.6307437E+02 2.00000006+00 C. 


ORACLS PRQGRAN 


4.O306688E+O2 9.848526XE+00 2.05XX054E+0X 6.6485X22E-01 


4.23564X4E+02 9.8496065E+00 2.1522331E+01 6.8264483E-0X 


4 .43263576+02 9. 8506871E+00 2.25337196+01 7.0044040E-0X 


4.6296516E+02 9. 8517678E+00 2. 35452X8E+01 7.1823790E-0X 



EXAMPLE 3 - MODEL-FOLLOWING 


QRACLS PRQGK^N 


2.66A1925S+00 -8. 4518962 fc-OX 


•5000000E+02 


7.8431432t+00 2. 5759893E +01 6.79Q0715E-01 4.e266891E402 9, 8526466E+00 2«4556828E^01 r«3603736E-01 

4.8277899E+02 2 .00000006+00 0. 

2.7042985E+00 -8.39101096-01 



Example 4 - Kalman-Bucy Filter 

This problem illustrates the asymptotic optimal estimator design capability 
of ORACLS. Consider the linear system 

x(t) = A x(t) + B u(t) + C(t) 


with A and B from example 1 and C(t) a zero-mean Gaussian white-noise 
process with constant intensity Q = Q’ ^ 0. If the system is digitally con- 
trolled using a zero-order hold with sampling times equal to integral multiples 
A = 0.05 second, the state equation becomes 

xjji+DA^ = A x(iA) + B u(iA) + 5(iA) (i = 0,1,...) 


where 


A = e^^ 


B 



dT 


and ^(iA) is a zero-mean Guassian discrete white-noise process with variance 
matrix 


Q 


^A 

j e-^^Qe^' dT 

'-'o 


Observations are made at the sampling instants and are modeled by 


y(iA) = H x(iA) + D u(iA) + ri(iA) 


where 



1.0 

0.0 

0.0 

0.0 

H = 

0.0 

1.0 

0.0 

0.0 


_-0.0852 

-0.0421 

-2.24 

-0.0021_ 
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0.0 

0.0 

0.0 

0.0 

0.713 

-0.403 


and ri(iA) is a zero-mean Gaussian discrete white-noise process with variance 
R = R’ > 0. Assuming ^(iA) and Ti(iA) are mutually uncorrelated, find a 
Kalman-Bucy predictor to asymptotically estimate the state x(iA) from knowl- 
edge of the measurements y(iA) and control inputs u(iA) up to time (i-1)A. 

If the estimate of x is denoted by x, the predictor has the structure 
(ref. 4) 

x[(i+1)Aj = A x(iA) + B u(iA) + F[y(iA) - H x(iA) - D u(iA^ 

The filter gain F can be computed directly from ASYMFIL after ignoring the B 
and D matrices__in the digital plant and observation equations. As in exam- 
ple 2, A and B follow from EXPINT. The matrix 5 can be found from SAMPL 
using A’ in place of A. Finally, we take Q and R to be 4 x and 
3x3 identity matrices, respectively. 

The executive program and output data follow. 


184 



progra^i kbfil 


7^/7A OPT«l 


FTN 4,6+439 


77/07/27, 08.37#18 


1 


5 


10 


15 


20 


25 


30 


35 


40 


PROGRAH KBFIL( INPUT# OUTPUT#! APE5-INPUT#TAPE6-0UTPUT) 

r 

DIMENSION A(16)#B (8)#H(12)#Q(l6)#R<9),ATIL(X6)#BTIL(e)#QTIL(l6)#F( 
112)# P(16)»DUMMY(114) 

DIMENSION NA(2)#NB(2)#NH(2)#NQ(2)#NR{2)#NATIL(2)»NBTIL(2)#NQTIt(2) 
l#IOP(5)#NF(2)#NP( 2)#H0UM(2) 

LOGICAL IOENT#OISC# NEWT# STABLE# FNULL 

C 

c 

C INPUT HOILFRITH DATA FOR TITLE OF OUTPUT 
CALL ROTITL 

C 

C INPUT COEFFICIENT MATRICES AND NOISE INTENSITIES 
C FOR CONTINUOUS SYSTEM 

CALL RE AD (5# A# NA# B# NB# H# NH# Q# NQ# R # NR ) 

C 

C GENERATE COEFFICIENT MATRICES FOR DIGITAL SYSTEM 

ItjP(l)«0 
OELT* .05 
Nl«(NA(l)*+2) + 1 

CALL EXPINT( A#NA# ATIL#NATIL#DUMMY#NA#OELT>IOP#DUMMY<Nin 
CALL MULT(OUMMY#NA#B#NB#BTIL#NBTIL) 

CALL PRNT(ATIL#NATIL#4HATIL# 1) 

CALL PPNT(8TIL#NBTIL#4H6TIL#1) 

C 

C COMPUTE VARIANCE MATRIX FDR PROCESS NOISE OF DIGITAL SYSTEM 
CALL c0l)ATE(0#N0# OTIL#NQTIL) 

CALL TRANP( A#NA#OUMMY#NDUM) 

I0P(2)« 0 

CALL SAMPL<DUMMY#NDUM#B#NB#OTIL#NQTIL#R#NR#W#NW#OFLT#10P#DUMMY(N1) 
1 ) 

CALL PRNT(QTIL#N0TIL#4H0TU#1) 

r 

C SOLVE FOP DIGITAL FILTER GAIN 

lOENT • .TRUE. 

DISC « .TRUE. 

STABLE » .FALSE. 

FNIJLL « .TRUE. 

NEWT « .TRUE. 

ALPHA*. 5 
I0P(l)»l 
inP(2)*0 


00 

U1 



00 


POOGRA^' KRcil 7<»/7^ 0PT«1 


FTN <»,6+<r3<) 77/07/27. 08.37.18 


<»5 


50 


I0P(3)-0 
IOP<^)»0 
IOP( 5)«0 

CALL ASYMFlL(ATIL^NATlL,G,N6,H#NH,Qm,NQTIL#R#NR,F,NF,P,NP,IDENT> 
IDI SC >NFWT>ST ABLf^FNULLf ALPHA. lOP^DUMHY) 


STOP 

tND 


EXAMPLE 


kalmas-pucy filter 


A MATRIX ROWS 

•2. 6000000 E+00 2 . 5000000E-01 

‘7,5000000E-02 -2 tTOOOOOOE-Ol 

7.0OOOOOOE-O2 -9.900C0O0E-O1 
l.OOOOOOOf+00 7.8000000F-02 

B MATRIX ^ ROWS 

l#7000000c+01 7.0000000E+00 

8*20000006-01 -3 *20000006+00 

0* A.6000000E-02 

0 * 0 . 

H MATRIX 3 ROWS 

i.oooooooe+00 o. 

0* l.OOOOOOOE+OO 

*8*5200000E-02 -A .2100000E-02 


0 MATRIX 

l.OOOOOOOE+OO 

0* 

0* 

0* 

R MATRIX 

l.OOOOOOOE+OO 

0 * 

0* 


A ROWS 

0 . 

l.OOOOOOOE+OO 

0 . 

0 . 

3 ROWS 

0. 

l.OOOOOOOE+OO 

0 . 


A COLUMNS 
-3.8000000E+01 
A.AOOOOOOE+OO 
-2.3000000E-01 

0 . 

2 COLUMNS 


A COLUMNS 

0. 

0 * 

-2.2A00000E+00 
A COLUMNS 

0. 

0* 

l.OOOOOOOE+OO 

0. 

3 COLUMNS 

0. 

0. 

l.OOOOOOOE+OO 


ATIL MATRIX 
8.7A60216E-01 
-3.0699872E-03 
3.77A31376-03 
A.6820759E-02 


A ROWS 

5.6206553E-02 

9.811A690E-01 

-A.8707923E-02 

A.9177112E-03 


A COLUMNS 
-1.76AA035E+00 
2.1998505E-01 
9.7958A37E-01 
-A.AB09768E-02 


BTIL MATRIX A ROWS 2 COLUMNS 

7.9692392E-C1 3 . 223A606E-01 

3.9250218E-02 -1. 5895733E-01 

6. 180065 AE-OA 6. 86761 37E-03 
2.0A35169E-C2 7. 90 520 A5E-03 


OTIL MATRIX A ROWS A COLUMNS 


00 


ORACIS PROGRAM 


0* 

0. 

5.2000000E-02 


C« 


0* 

0 . 

-2.1000000E-03 


0 * 

0* 

0 . 

l.OOOOOOOE+OO 


-2.352A288E-03 

2.8616A63E-0A 

2*5772852E-03 

9.9996068E-01 



00 

00 


EXA'IPLF A KALMAN-BIXY FILTGP 


ORACIS 


9.706<»533E-O2 

-5*6349895f-03 

‘A*<r6^3769c-02 

2.0737669E-03 


-5*63^9698E-03 

<»#9959106t-02 

^•2255^961-03 

-1.21695det-05 


-4».<»6^3769E-02 

<>•22554966-03 

4*91721«4E-02 

-6^?173577E-04 


2.0737669E-03 

-1^2169539E-05 

-6.ei73577£-04 

5*0057551E-02 


PROGRAM TO SOLVE THE DISCRETE INFINITE-DURATION OPTIMAL FILTER PROBLEM 


A MATRIX 
8.7460216E-01 
-3^069B372E-03 
3^7743137E-03 
4.6820759E-02 


4 ROWS 

5^6206553F-02 

9,eil4690E-01 

-4.8707923E-02 

4.91771126-03 


4 COLUMNS 
-1^7644035E+00 
2.199e505E-01 
9.79584376-01 
-4.48097666-02 


-2.3524208E-O3 

2.8616469E-04 

2.5772852E-03 

9.9996068E-01 


G IS AN IDENTITY MATRIX 


H MATRI 

i.oooooooe +00 

0 • 

-6.52000006-02 
INTENSITY MATRIX 


X 3 ROWS 

0. 

X.OOOOOOOE+00 

-4.2100000E-02 

FOR COVARIANCE 


4 COLUMNS 

0. 

0 . 

-2. 24000006+00 
OF MEASURtMENT 


0« 

0 . 

-2.1000000E-03 

USE 


R matrix 3 ROWS 

i.oooooooe +00 o. 

0. l.OOOOOOOE+00 

0 . 0 . 

INTENSITY MATRIX FOR COVARIANCE 


Q MATRIX 
9.7064533E-O2 
-5.6349895E-03 
-4.4643769E-02 
2.07376696-03 


4 ROWS 

-5.6349895F-03 

4.9959106E-02 

4.2255496E-03 

-1.2169588E-05 


3 COLUMNS 

0 . 

0 . 

l.OOOOOOOE+OO 
OF PROCESS NOISE 


4 COLUMNS 
-4.4643769E-02 
4.2255496E-03 
4.9172184E-02 
-6.8173577E-04 


2.0737669E-03 

-1.2169588E-05 

-6.6173577E-04 

5.0057551E-02 


PROGRAM 


FILTER GAIN 



tX^MPtE A 


kalmam-hucy filter 


OPACLS PROGRAM 


F *<ATRIX A ROWS 3 COLUMNS 

5.3^5^2336-01 -1.00008 70t-02 3. 6 31 271 7E-01 

-^.36530816-02 1 • 8601P 1 ^c-01 -3 . ^r32 ^ A6 « E-02 

-8.A771306E-02 -9.0256^6^6-03 -1.21357066-01 

-7.O9270<f2E-O3 3 . 8 1 65991 F-02 -1 . 9 1 1 375 9E-01 

STEADY-STATE VARIANCE NATRIX CF R ECCINSTPUC T I ON ERROR 


P MATRIX ROWS ^ COLUMNS 

8.86496806-01 -5.71943366-02 -2 . 2 09402 5 c-01 -1 . 26 260 90E-01 

-5.71943866-02 2.35575736-01 8.303025^(6-03 5. 2608810E-02 

-2. 209402 5E-01 8 . 30302 5R E-03 1. 056 5196E-01 9. 1725042E-02 

-1.262 6090 E-01 5.28086106-02 9. 1 72 5C4 2E-02 1 . 9945 24 8E+01 

EIGENVALUES OF P 


EVLP MATRIX 4 ROWS 1 COLUMNS 

4.72574956-02 
2.3111030E-01 
9.4813881E-01 
1.99466666+01 

A-FH MATRIX 4 ROWS 4 COLUMNS 

3.70993276-01 6. 14953 86E-02 -9. 5099664 F-01 -1.58986176-03 

3.76587496-02 7. 93603 656-01 1 . 4 30962 4 E-01 2. 1408330E-04 

7.02059986-02 -4. 4791409E-02 7. 0 774 4 5 5 E-01 2 . 3 2 24 3 54E-03 

4.18866216-02 -3. 91 901 72 E-02 -3. 609579 7E-01 9.99664296-01 

EIGENVALUES OF A-FH MATRIX 


5. 40242 396-01 ? . 2 217840E-01 

5.4024239E-01 -2. 221 76 40F-01 

7.9425006E-01 0. 

9.97355926-01 0. 


00 

VO 


1 



CONCLUDING REMARKS 


This report has presented some 60 subroutines (43 primary purpose and 
17 supporting) which can be used for the analysis and design of state-variable 
feedback control laws for time-invariant linear systems. The basic synthesis 
approach is through the well-founded linear-quadratic-Gaussian (LQG) method- 
ology. The ORACLS system represents an attempt on the author's part to employ 
some of today's best numerical linear algebra techniques to the LQG design 
problem. A modular structure has been Incorporated so that new methods can be 
incorporated, as they are developed. 

The examples presented indicate only some of the capability of the ORACLS 
package. The program has been applied at the Langley Research Center to prob- 
lems in aircraft design using state- variable equations of order up to 16. The 
Liapunov equations (subroutines BARSTW and BILIN) have been used successfully 
to solve equations of 79th order. 

Each subroutine was compiled on a FORTRAN Extended compiler and executed 
on a Control Data Cyber digital computer system. The full ORACLS package 
(60 subroutines) requires 22 345 decimal (53 511 octal) words of storage. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
November 1 , 1977 
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